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LIS INTRODUCTION 


A. BACKGROUND 

The field of cyclic spectral analysis is growing in 
importance as non-stationary modulation schemes proliferate in 
modern communications. Traditionally, spectral analysis 
methods have assumed that signals of interest were stationary 
or had  non-time-varying statistical parameters. This 
assumption is untrue for the cyclostationary signals created 
when using pulse or carrier amplitude modulation, phase or 
frequency carrier modulation or digital pulse or carrier 
modulation. Gardner [Ref. l:pp 419-453] discusses many of 
these signal types. Spread spectrum modulation techniques 
such as direct sequence PSK and frequency-hopped FSK also 
exhibit temporally varying spectral features [Ref. l:pp 453- 
457]. Reference 2 provides a tutorial on cyclic spectral 
analysis. 

Cyclic spectral analysis methods take advantage of 
spectral variation over time by correlating spectral 
estimations at discrete time intervals to locate signals that 
may not otherwise be identified. A direct sequence signal 
with a very wide bandwidth may be indistinguishable from a 
noisy background in the frequency spectral estimate of 


Figure 1 (Réf. 3:pp 2-9 - 2-10]. But this signal becomes 


clearly identifiable in the cyclic spectral estimate of 
Figure 2, 

One cyclic spectral technique is the frequency smoothed 
cyclic periodogram method (FSM) [Ref. 1:pg 464]. This has 
been implemented in the Cyclic Spectral Analysis Software 
Package [Ref. 4] as the program sxaf. Hereafter, this 
software will be referred to as the SSPI package or the FSM. 
This method was used to generate Figures 1 and 2. 
Unfortunately, the FSM is less efficient than time smoothing 


methods for general spectral estimation [Ref. 5:pg 38]. 


B. THESIS GOALS 

The purpose of this thesis is to implement two time 
smoothing algorithms, the FFT Accumulation Method (FAM) 
[Ref. 5:pp 44-47] and the Strip Spectral Correlation Algorithm 
(SSCA) [Ref. 5:pp 47-48]. The FAM is less computationally 
complex than the FSM and the SSCA less so than the FAM. Both 
methods have been written to be compatible with the SSPI 
package. By integrating FAM and SSCA into the SSPI package 
future research may more easily be done in algorithm 
performance comparison and enhancement. 

Both the FAM and SSCA accept input data generated by the 
SSPI utility progran bau They also generate results 
compatible with the SSPI plotting utility plot sarita 
program may optionally generate output in a format compatible 


with the MATLAB [Ref. 6] package which is in common use at the 


Naval Postgraduate School. Input and output data are in ASCII 
file formats. 

Throughout this thesis the same BPSK signal data is used 
to illustrate the implementations of FAM and SSCA. A 
representative signal sample is plotted in Figure 3. No noise 
Was added to this signal. By estimating the frequency 
Spectrum at different points in time, it is evident in Figure 
4 that the spectral features do indeed vary temporally. 

Chapter II describes the implementation of the FAM. 
Chapter III describes the implementation of the SSCA. Chapter 
IV describes a utility to estimate cycle frequencies for a 
given baseband frequency. The Sub-FFT Accumulation Method 
(SUBFAM) program [Ref. 7] gives a quick result for a specific 
subset of the spectral plane. Chapter V discusses specific 
algorithm performance comparisons. Finally, Chapter VI 
presents conclusions and recommendations for future areas of 


research. 
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Figure 1 Power Spectral Density of BPSK Signal in Noise 


(Ref. 4:pg 18] 





Figure 2 Cyclic Spectral Estimate of Direct Sequence Signal 


in Noise [Ref. 4:pq 18] 


P 





Figure 4 Time Sequence Spectrum of Typical BPSK Signal 
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II. FFT ACCUMULATION METHOD 


A. THEORY 

The FFT Accumulation Method (FAM) [Ref. 5:pp 44-47] is a 
Fourier transform of correlation products between spectral 
components smoothed over time. Periodicities in the spectral 
components then become detectable. The cyclic spectral plane 
as shown in Figure 5 ranges in frequency from -f,/2 to +f,/2, 
and in cycle frequency from -f, to «f, where f, is the sampling 
frequency. For each unique(f,&) cermin krg Be Er 


cross-spectral estimate is defined to be [Ref. 5:pg 44]: 





P-1 -i2nrq 
+gqA : = 
n (on ee X (ELLE) Y? (EL, E nm (1) 
r=0 


The FAM auto-spectral estimate is defined to be 


[Ret . 5:00 241 


U 


-1 -i2mrga 
a,+gha DEN LE 


Six, (nL, f,)4, X.(rL f.) X; ie f (2) 


H 
II 
° 


where: j is the average frequency bin (k+f)/2 
q is the difference or bandwidth frequency (k-/) 


g.(n-r) represents the Hamming window. 


The FAM was implemented by forming an array from x(kT) 
SS Nel) with rows which are N points long from the 
input sample data. The starting point of each succeeding row 
is offset from the previous rows starting position by L 
samples. A Hamming window is applied across each row which is 
then Fast Fourier transformed and downconverted to baseband. 
MM ul GE this point is a two-dimensional array with 
columns representing constant frequencies. Each column was 
point-wise multiplied in turn with the conjugate of the 
columns resulting from processing y(kT) (0 s k s N-1). Each 
resultant product vector was Fast Fourier transformed and the 
low frequency half placed into the final cyclic spectral plane 
at the appropriate location in the cyclic spectral plane. 
Figure 6 shows a block diagram for the FAM cross-spectral 
estimate. 

The following subsections in this chapter discuss in 
detail what has been described in the previous paragraph. The 
cycle frequency resolution of FAM is Aw = f,/PL [Ref. 5:pg 44] 
where f, is the original data sampling rate. The frequency 
2 Ol Ono FAM is Af = £/N‘ [Ref. 5:pg 44]. The time- 


sI ie rc Sol ton product is AtAf - PL/N' [Ref. 5:pg 45]. 


The command line format for calling the FAM program is 
provided in Appendix A. The FAM source code listing is 


provided in Appendix B. 


B. IMPLEMENTATION 
1. Input Channelization 
The input sample data is formed into a two-dimensional 


array. The array row length is equal to the number of input 


Channels N’. For a given number of input sample points N, a 
row size of N', and a chosen offset L, there are P = (N- 
N')/L = N/L rows formed. The choice of N' must take into 


consideration that ideally the time-frequency resolution 
product must be much greater than one [Ref. 5:pg 40]. N' 
should also be a power of 2 to avoid truncation or r 
padding in the FFT routines. L should be chosen to be less 
than or equal to N'/4. 

The completely filled array is P rows by N' columns. 
Figure 7 shows how a small array is filled from a discretely 
sampled signal x(kT) when N'=16, P=8 and L=4. The number 
inside each cell represents the value of k used to index on 
x(KT) to fill that location in theWarrad 

Figure 8 shows the magnitude of the original data for 
the example BPSK signal organized into the P by N' array where 
N=512, N'=32, L=4 and P=128. The input data is assumed to be 


complex with a real and imaginary component to each sample 


point. Figure 9 shows a single row of the array. The phase 
changes of the BPSK are evident. 
2.  Windowing 
A Hamming window [Ref 7:pg 467] is applied to each row 
of the array. The equation for the Hamming window is: 


ent, 


EU a (3) 
nus 


wea 540 d6cost 





A 32-point Hamming window is plotted in Figure 10. It 
is applied to both the real and imaginary parts of the complex 
example array. The magnitude of the resultant array is shown 
in Figures 

3. First FFT 

Each row of the windowed data array is Fast Fourier 
transformed to reveal the first spectral components. The 
resultant array is still indexed P rows by N' columns but now 
the column index relates to a specific bin of spectral 
frequencies. Figure 12 illustrates this relationship. 
Figure 13 shows the results of FFTing the BPSK example. 

4. Downconversion 

Each row of spectral components is downconverted to 

baseband through multiplication with the complex exponential, 


-i2nkmL 
e 


where: m is the row index, O0 sm s P-1 


k is the column andex, Os "a< 


The magnitude of the exponential is unity over the 
array but the phase shows considerable variation. Figure 14 
shows the phase of the exponential over the P by N array. 
Figure 15 shows the phase for one representative row. The 
magnitude of the array remains unchanged from Figure 13. 

5. Multiplication 

For each cell in the cyclic spectral plane, one column 
of the array is multiplied with the conjugate of another 
column. The separation of the columns is determined by the 
desired frequency separation or cycle frequency (a=f,-f,). The 
mid-point between the columns is the frequency bin of interest 
(f= (ot oe Figure 16 shows two columns to be conjugate 
multiplied for use in filling a specifie f,@ eellz E1 o 
illustrates the f,&œ bin values in the cyclic spectral plane 
box =< | Figure 18 shows the corresponding two column 
indices used to form the product vector which is used to 
produce the cells of Figure 17. 

6. Second FFT 

The product from the previous multiplicatilen is FRI d 
to yield a P-point result. Only the middle of the resulting 
spectrum is retained and stored into the designated f,a cell. 


The upper and lower ends are undesireable because of increased 


10 


estimate variation at the channel pair ends [Ref. 5:pg 45]. 
Figure 19 illustrates which part of the estimate is retained 
and placed into the cell. Figure 20 shows the final result of 
the FAM auto-spectral process on the example BPSK signal. 

7. Data Reduction 

The FAM program typically generates large output data 
files. For convenience, an option may be chosen to reduce the 
Eneumt of outpub: By comparison sorting for the largest qo 
value in an f,o« cell, the number of a slices output is reduced 
Pon N'B/2 “to. Nî Overall FAM program execution time 
increases accordingly to accomplish the sorts. 

Figure 21 illustrates the output full cyclic spectral 
plane storage array before sorting. Figure 22 shows the 
output array after data reduction has been completed. Figures 
23 and 24 plot the resulting spectral half-planes without and 
with data reduction respectively. 

Since all cells have an equal number of a values, all 
sorts are of equal complexity. Further work on data reduction 
would require selection of a largest o value from widely 
varying numbers of «a values depending on the choice of N, N' 


and L. 


C. PERFORMANCE 
An established method of evaluating the complexity of an 


algorithm is to determine the number of floating point 


qe 


operations that must be performed. Fox this estimate it is 
assumed that an auto-spectral estimate is being computed and 
the data is complex in every step. It is also assumed that 
the Hamming window coefficients and the downconversion 
coefficients have been previously computed and stored for 
later use. Each N-point FFT requires (N/2)*log,N complex 
floating point multiplies [Ref. 8:pg 506] or 4*(N/2) *log,N 
real floating point multiplies. The costiogt any mouc pie anna 


reduction is not considered here. 


Applying the window ee. —— 2*P*N' 

First FET. IO AA NN 2*P*N' *log,N' 
DownconversXono-.-. See ee E E E EA, 4*P*N' 
Múltiplications HS. ANN 4*N'*N'*P 

Second FEAR 2*N'*N'*P*]og,P 
Total:  (6«4*N')*P*N' r (2*P*N')*(log,N'  N'*log,P) (4) 
Since P=N/L (5) 
Total:  (6«4*N')*NN'/L + (2*NN"/L) (1log,N” + N'log,N/L) (6) 


12 
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Figure 6 Block Diagram of FAM Cross-spectral Estimate 
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Figure 7 Layout of A Sample Input Array Showing Input 
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Figure 16 Generic Array Showing Two Columns to be 


Multiplied, f = -1, œ = 2 
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Generic Array Showing Column Index Values 
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Figure 19 Retained Section of the Second FFT Results 
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Figure 20 Final FAM Results of the BPSK Example, 


N=512, N'=32 and L=4. 
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FAM Storage Before Data Reduction 


Figure 21 
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| Figure 22 FAM Storage After Data Reduction 


25 






= 0 5 


À 


1.92e+03 


Figure 23 FAM Result Without Data Reduction, 


N=512, N'=32 and L=4. 
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Figure 24 FAM Result With Data Reduction, 


N=512, N'=32 and L=4 
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III. STRIP SPECTRAL CORRELATION ALGORITHM 


A. THEORY 


The Strip Spectral Correlation Algcorsa hm (SSCA) 
[Ref. 5:pp 47-48] is a Fourier transform of eorrelacrıen 
products between spectral and temporal components smoothed 
over time. Periodicities in the spectral components then 
become detectable. The cyclic spectral plane as shown in 
Figure 25 ranges in frequency from -f,/2 to +f,/2, and in cycle 
frequency from -f, to +f, where f, is the sampling frequency. 
For each unique (f o mse io pan Figure 25, Ehe SSCAMCLOSON 


spectral estimate is defined to be [Ref. 5:pg 47]: 


f,*gâa 


£ -i2xrq 
Sxy. (n I gAa 


mar 
de 2 M EISE r ac B G) 





The SSCA auto-spectral estimate is defined to be 


(Ref. S:pg 47]: 





f,*gâa fa qå a Pei E 
xx, A Foa- ae, Ac Ur a n O j (8) 
r=0 
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where: k is a multiple of the frequency resolution, 
“Ney 2 NIE TONNES 
q is a multiple of the cycle frequency resolution, 
WE GW NF) 


g(n-r) represents the Hamming window. 


The SSCA was implemented by forming an array from x(kT) 
MES E < N-2) with rows which are N’ points long from the 
input sample data. The starting point of each succeeding row 
is offset from the previous rows starting position by L 
samples. A Hamming window is applied across each row which is 
EnenaewastoWourier transformed. and downconverted to baseband. 
The result at this point is a two-dimensional array with 
columns representing constant frequencies. Each row is 
transposed and replicated L times for a total of PL columns. 
Each column is then multiplied by a sample y(kT) (0 < k s (PL- 
D) Each resultant row of the array is Fast Fourier 
tonstorwmed and placed into a strip of the final cyclic 
spectral plane at the appropriate location. Figure 26 shows 
a block diagram for the SSCA cross-spectral estimate. 

The following subsections in this chapter discuss in 
detail what has been described in the previous paragraph. The 
cycle frequency resolution of SSCA is Aw = f,/N [Ref. 5:pg 48] 
where f, is the original data sampling rate and N is the 


number of input samples used. The frequency resolution of 
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SSCA is equal to the sampling rate divided by the number of 
channels N’, Af = £,/N' [Ref. 5:pp 42,48]. The timestrequeney 
resolution product is AtAf = N/N' [Ref. 5:pg 48]. 

The command line format for calling the SSCA program is 
provided in Appendix C. The SSCA source code listing is 


provided in Appendix D. 


B. IMPLEMENTATION 
l. Input Channelization 
The input sample data is formed into a two-dimensional 


array. The array row length is equal to the number of input 


channels N'. For a given number of input sample points N, a 
row size of N’, and a chosen offset L, there are P = (N- 
N')/L = N/L rows formed. The choice of N' must take into 


consideration that ideally the time-frequency resolution 
product must be much greater than one [Ref. 5:pg 40]. N' 
should also be a power of 2 to avoid truncation or zero- 
padding in the FFT routines. L should be chosen to be less 
than or equal to N'/4. 

The completely filled array is P rows by N' columns. 
Figure 27 shows how a small array is filled from a discretely 
sampled signal x(kT) when N'=16, P=8 and L=4. The number 
inside each cell represents the value of k used to index on 
X(kT) to fill that location in the array. 

Figure 28 shows the magnitude of the original data for 


the example BPSK signal organized into the P by N’ array where 
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N=512, N’=32, L=4 and P-128. The input data is assumed to be 
complex with a real and imaginary component to each sample 
point. Figure 29 shows a single row of the array. The phase 
Changes of the BPSK are evident. 
2. Windowing 
A Hamming window [Ref 7:pg 467] is applied to each row 


of the array. The equation for the Hamming window is: 


E Du Me 1 (9) 


AR SA .46C0Ss 4 





A 32-point Hamming window is plotted in Figure 30. It 
is applied to both the real and imaginary parts of the complex 
example array. The magnitude of the resultant array is shown 
Un Elqure 3c 

s PIrsSt FET 

Each row of the windowed data array is Fast Fourier 
transformed to reveal the first spectral components. The 
resultant array is still indexed P rows by N' columns but now 
the column index relates to a specific bin of spectral 
frequencies. Figure 32 illustrates this relationship. 
Figure 33 shows the results of FFTing the BPSK example. 

4. Downconversion 
Each row of spectral components is downconverted to 


baseband through multiplication with the complex exponential, 


Spl 


where: m is the row index, O s m s P-1 


k is the column zndex, O WINNIE 


The magnitude of the exponential is unity over the 
array but the phase shows considerable variation. Figure 34 
shows the phase of the exponential over the P by N' array. 
Figure 35 shows the phase for one representative row. The 
magnitude of the array remains unchanged from Figure 33. 
5. Replication 
Each row is copied into one column of an empty N' Dy 
PL array. It is then replicated in the L-1 adjacent columns. 
Figure 36 illustrates this process. 
6. Multiplication 
Each column of the array is pointwise multiplied with 
the conjugate of a sample value y(kT). There are PL=N columns 
and PL samples from y(kT). Figure 37 illustrates the 
conjugate multiplication process. 
7. Second FFT 
Each row from the previous multiplication is Fast 
Fourier transformed to yield a PL-point result. Each 


resulting vector is stored into a strip of the cyclic spectral 


Br 


plane. Figure 38 shows the final result of the SSCA auto- 
spectral process on the example BPSK signal. 
8. Data Reduction 

The SSCA Program typically generates large output data 
files. For convenience, an option may be chosen to reduce the 
amount of output. By comparison sorting for the largest a 
value in an £f,& cell, the number of a slices is reduced from 
N to N/L. Overall SSCA execution time increases accordingly 
to accomplish the searches. 

Figure 39 illustrates the output full cyclic spectral 
plane storage array before sorting. Figure 40 shows the 
output array after data reduction has been completed. Figures 
41 and 42 plot the resulting spectral half-planes without and 
with data reduction respectively. 

Since all cells have an equal number of a values, all 
sorts are of equal complexity. Further work on data reduction 
would require selection of a largest o value from widely 
varying numbers of a values depending on the choice of N, N' 


and L. 


C. PERFORMANCE 

An established method of evaluating the complexity of an 
algorithm is to determine the number of floating point 
operations that must be performed. For this estimate it is 
assumed that an auto-spectral estimate is being computed and 


the data is complex in every step. It is also assumed that 
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the Hamming window coefficients and the downconversion 
coefficients have been previously computed and stored for 
later use. Fach N-point FFT requires (N/2)*log,N complex 
Floating point multiplies Ref. 8:p9g 506P Or 2*N*log N Ee 
floating point "multiples The cost of any output dara 


reduction is not considered here. 


Applying the windows ` sp AN 2*P*N' 

First FET. s ss 2*P*N'*log,N' 
Downconversien.. Mi...  — A*P*N' 
Multiplication V —— NEUE 4*N*N’ 

Second EBTI Ae i aeaea ee ee 2*N*N'*log,N 
Total: 2*N’*(6*P + 4*N + (2*P + 2*N) *log,N) 61:0) 
since P=N/L (119 
Total: 2*N’*((6*N/L + 4*N) + (2*N/L + 2*N) *log,N) (12) 
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Figure 25 Generic SSCA Cyclic Spectral Plane 
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Figure 26 Block Diagram of SSCA Cross-Spectral Estimate 
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Figure 30 Thirty-two Point Hamming Window 
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Figure 37 Multiplication of Columns by y* 


44 






A 
WAI e 
D U Ly ES ì 

HA qd Me 

| SR “ Al ff It, | Ke, A 
| S 941 = y) ; IA JUN ŷ IN; I 

4 77 ud A WY 

"I, J 

ARISTAS 


iub e ud e /// 1/41 t 
NIZA ADE ANA) // A | 
TIN GDR o A | 
I DAE I EE GWA 


-0.5 


0.5 
479. 


Figure 38 Final SSCA Results of the BPSK Example, 


N=512, N'=32 and L=4. 
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T : 
j Storage Array Before Data Reduc 
Figure 39 SSCA © NOR 
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Figure 40 SSCA Storage Array After Data Reduction 
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Figure 41 BPSK Example Results Without Data Reduction, 


N=512, N'=32 and L=4. 
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IV. SUB-FFT ACCUMULATION METHOD 


A. THEORY 
In the course of this thesis work a need arose toma 
flexible and customized Spectral Correlation Function (SCF) 


[Ref. 7] to determine the presence of phase-shift keyed 


Signals in experimental data. The SCF  [Ref. 7] is 
represented: 
N-1 an (Erz, 
Saale AOS eee | (13) 
r=0 
where: 


N is the number of data samples used 
f; is the intermediate modulation frequency 
f, is the sampling frequency 


OSF is a one-sided bandpass filter 


The program written to implement this algorithm is named 
the Sub-FFT Accumulation Method (SUBFAM) program to 
distinguish it from the more general purpose FAM program. The 
command line format for calling the SUBFAM program is provided 
in Appendix E. The SUBFAM source code is provided in 


Appendix F. 
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B. IMPLEMENTATION 

Figure 43 illustrates the operations performed in the 
SUBFAM program. A set of N points is obtained from within a 
file of signal data values. Each set of N points is processed 
through a one-sided bandpass filter to remove either all 
positive or negative frequencies as desired.  Downconversion 
of each file to baseband from the original sampling and 
transmission frequency bands follows. The results are then 
correlated through a complex conjugate multiplication. A 
final N-point Fast Fourier Transform yields the spectral 


correlation result. 


C. PERFORMANCE 

The SUBFAM program performed as required. Figure 44 
illustrates the results of processing test data containing 
phase-shift-keyed signal samples. The peaks at the chip 
frequencies are apparent. Figure 45 shows the results of 
correlating test data with data containing only white noise. 
The lack of correlation between spectral features yields 
results which are four orders of magnitude less than in the 


successful signal detection of Figure 44. 
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Figure 43  SUBFAM Program Block Diagram 
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V. ALGORITHM PERFORMANCE 

The complexity of the FAM and SSCA algorithms is estimated 
MERD per II Section C and Chapter III Section C 
respectively. Both require on order of Nlog.N floating point 
multiplies where N is the original number of sample points 
used. Brown and Loomis [Ref. 9] explore the complexity of the 
FAM, SSCA and FSM algorithms in some detail. Their results 
are consistent with an order Nlog,N multiplies for FAM and 
SSCA. 

The cyclic spectral plane is symmetric about the f=0 and 
the a=0 axes for real signals. For the most efficient 
performance it is only necessary to compute one quadrant of 
the cyclic spectral plane. To compare the performance of FAM 
and SSCA with FSM, complexity figures derived from Reference 
9 for quadrant complexity will be used. The three algorithms 


require the following numbers of floating point multiplies: 


Pai Complexity (Ref. 9:pg 18] 


2NN'10g, (57) + 8Nlog, (N') + 4NN' - 20N (14) 
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SSCA Complexity Ret. 9:pg 20) 


NN’log,(N) + 2NN’ + 8Nlog,(N’) + 12N (359 


FSM Complexity |Re£. 9 pg 16] 


N* + ANI GG (NN) (16) 


The complexity of four example cases with varying N and N' 


are estimated below using eguation (14), (15) and (16). 
N N' FAM SSCA FSM 
1024 64 5.3*10° 3.6*10> 105 

2048 128 3 0210 AOS 4 2x10 
32,768 512 158109, 7.1319 20, 
1,048,576 8192 8.0*10!9 TO LON 10 


It is evident that FAM and SSCA require substantially 
fewer multiplies than FSM. Particularly with larger values of 
N, FAM and SSCA require Nlog,N while FSM requires N° floating 


point multiplies. 
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VI. CONCLUSIONS 


A. SUMMARY 

The purpose of this thesis was to implement the FAM and 
SSCA for use in cyclic spectral analysis researcn. Both 
algorithms were successfully implemented in a manner 
compatible with the SSPI analysis package to facilitate future 
work. By allowing the user a choice of two output file 
formats, results may easily be used in either MATLAB [Ref. 6] 
or SSPI [Ref. 4] for further signal processing applications. 

It has been shown that the FAM and SSCA algorithms 
generate results consistent with FSM but require substantially 
fewer floating point multiplies than FSM. This is especially 


true as the number of sample data points used increases. 


B. RECOMMENDATIONS 

In the course of this thesis it has become apparent that 
the FAM and SSCA algorithms consist of inherently parallel 
operations. Roberts and Loomis explore these possibilities in 
Reference 10. The high-level sequential nature of these 
algorithms also lend them to pipelining architectures. By 
combining a parallel or multi-processor approach with 
pipelining a large improvement in performance should be 


obtainable. This would be espectially true in applications 


5 


requiring large sets of signal data samples or repeated, 
time-sequenced cyclic spectral plane estimations. 

Two areas which could benefit from further work are 
automatic signal identification and signal parameterization 
algorithms. The results from cyclic spectral analysis 
algorithms such as FAM and SSCA may be input into Linear 
Associator or Mapping neurocomputer networks as described in 
Reference 11. Intensity transformations using digital image 
processing techniques from Reference 12 also appear to have 
promise. Both classes of techniques could have utility in 
automatic signal identification and parameterization. They 
also have the added benefit that they lend themselves to 


parallel and pipelined computing architectures. 
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Correct 


for the 


fam 
fam 


or for 


APPENDIX A. FAM PROGRAM USE 


commands are: 
cross-fam 


pure lino tale outputfile n nprime l Oflag 
inputfilel inputfile2 outputfile n nprime 1 Oflag red 


the auto-fam 


fam inputfilei Outputfile n nprime 1 Oflag 
gameonpütfileizaurpmerile n nprime L Oflag red 


where: 


more: 


input f 


inputfilel is the file containing the x signal samples 
inputfile2 is the file containing the y signal samples 
outputfile is where the spectrum values will be placed 
n is the number of samples to use from the inputfiles 
and it must be a power of 2 
nprime is the group size of the input datasets and it 
must be a power of 2 
1 is the starting offset of subsequent datasets and 
it must be a power of 2 
Oflag indicates the output filetype, ascii or binary 
-asc for an ascii, MATLAB compatible file, 
full cyclic spectral plane 
Po fora plot sxa, SSPI compatible file, 
half cyclic spectral plane 
red reduces the amount of data output 


the first two entries in every input file are expected 
to be the datatype and n. datatype = 1 for real 
values, 2 for complex values. n is the number of 
samples contained in the file, one sample per line. 
ile format: 

datatype n 


sample #1 
sample #2 


sample #n 
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output “file “ASCII MATLAB formar: 
nprime n 


OUtPpuUE E Lm 
output (LIZ 


output (nprime) (n) 


output file ' SSPI plOc SS armar. 
datatype (l For real) 2 ftorocompiex) 
number of alphasalpha minalpha max 
number of freqsfreq minfreg max 
value of alpha #1 


number ot tregs Jn WPEPedgem rcc 
spectrum at freq min 


Spectrum at ired max 


value of alpha #alpha_max 
number of freqs in #alpha maxfreq minfreq max 
spectrum at freq min 


spectrum at freq max 
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APPENDIX B. FAM PROGRAM LISTING 
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Dee 


9 12 2710819928 ram, fam espaces 


include «stdlib.h» 

zinclude «stdio.h» 

zinclude «math.h» 

sinclude "/home3/carter/thesis/SSPI/pam/fft.c" 
sinclude "/home3/carter/thesis/SSPI/pam/radix.c" 
main (argc,argv) 


T) 


argc; 


Char cargui]; 


{ 


COMPLEX. "Xx, *s3, y3, Gwneen Md 
COMPLEX *tempf£ft,*yl,*Sl,**S2, y2, Wu Ow 


/ * 
/ 


int 
Int 
ine 
int 


x passes data to/from fft routine 
ss holds results OP EEIuUUDOQESS 
y3 holds results cl Fr ling S 
dwnconv holds downconversion coefficients 
s4 holds correlation multipiyresunies 
tempfft passes data to/from fft routine 
yl holds y samples from inpuerile 2 
sl holds x samples from inputfile 1 
s2 holds channelized x samples 
y2 holds channelized y Samples 
window holds Hamming window values 
ER 
i,j,k,type,n,file n,nprime,p,l,curectiono norm; 
a,b,c,cross,num f,max num f,data type,max num alf; 
tempi, temp], halip, reduce, cure art, 
*outint; 


float numreal,numimag,convfac,mainfac,num alf,f min,f max; 
float alpha min,alpha max,f min all,f max all; 

float tempreal,tempimag,bigmag,tempmag; 

float twopai-6.2831858530705; 

float *outreal; 

double z, y; 


char 


*intilel, *infilo2, *outfrile; 


FILE rp, ::p2 SED: 
/* check for the correct number of input arguments 


the correct commands are: 


fam inputfilel inputfile2 outputfile n nprime 1 Oflag reduce 
fam ET inputfile2 outputfile n nprzume IC aq 

fam eee outputfile n nprime 1 Oflag reduce 

fam E Sutputfile n npeıme I lau 


where:  inputfilel is where the x signal samples are expected to be 


inputfile2 is where the y signal samples are expected to be 
outputfile is where the specttum values will be pute 
n is the number of samples to use from the inputfiles 
nprime is the group size of input datasets 
l is the starting offset of subsequent datasets 
Oflag indicates the output filetype, ascii or binary 
-asc for an ascii, matlab compatible file, full plane 
-plo for a plot sxaf compatible file; halt (plane 
reduce is an option to reduce the amount of output 


be 


eee 12:21 1992  fam/fam.c Page 2 


EE ^ Dco 3nd pn re expected o ce at the top otf che input file 
datatypen OER real,.,2 tor complex 
pnaxsscneenampersof samples following 

Beten by LEDR Nancy J. Carter, USN NEGS Monterey CA Ol October 1992 


powUnolement the Frequency Accumulation Method of cyclic spectral analysis 
as described in the Brown/Gardner/Loomis IEEE paper 


ANDE ctober 1992.... changed to asciz MATLAB output format 
EN ber 1992....added plot _sxaf output format compatible with SSPI 
EN vmber 1992... .added option to reduce amount of output 
Ez 
/* CHECK INPUT VALUES E 
Jut Es 
if ((argc !- 7)&&(argc != 8)&&(argc !- 9)){ 
pentf("fam....fatal errorMn"); 
EuSER("....... incorrect number of calling arguments\n"); 
E co correct formats are:\n"); 
ees | o ae um O purerfmle2zeutputfile n nprıme 1 Oflag reduce\n"); 
ALE acc o fam inputfilel inputfile2 outputfile n nprime i Oflag\n"); 
NEL (coc £am input ile outputfile nn nprime 1 Oflag reduce n"); 
FE E(",........ fam inputfile outputfile n nprime 1 Oflag\n"); 
EXSER(U.......... where\n"); 
BER (".......::: inputfilel contaıns the x samples\n"); 
Sete (wa tee inputfile2 contains the y samples\n"); 
Bentert(".......... eutputfile will contain the resultsin"); 
coo fois Chesnumoer or input isamples treo uSeAn"); 
EE... o... . nprime is the group size of input datasets Mn"); 
Euro . . ..... ... . l is the offset of subsequent datasets\n"); 
DE ccoo CF GIS the cuetpue £ile format \n"); 
e Oflag = Gasc, an ascii file 1s produced\n”),; 
FE”... ....... Grlag -=plo, orplotzsszar Elle Ts producedxn$j5 
Crab PNE reduce is an option for reduced amount of outputin"); 
EXTE(); 


) 


EE rgc--7) ( /* this iS an auto-fam, no output reduction */ 


cross=0; 
infilel=argv[1]; 
outfile=argv[2]; 
n=atoi (argv[3]); 
nprime-atoi(argv[4]); 
l=atoi(argv[5]); 
reduce=0; 


} 


if ((argc==8)&&(argv[7][0]!='r”')) 1 »/* this is a cross-fam, no output reduction  */ 


cross=1; 
infilel=argv[1]; 
infile2=argv[2]; 
outfile-argv[3]; 
n=atoi (argv[4]); 
nprime=atoi (argv[5]); 
l=atoi(argv[6]); 
reduce=0; 


} 


“rc «S(argvi?l[0]—— r'/)) 1 /* this ism auto-fam, with output reduction */ 


cross=0; 
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infilel-argv “MW: 
outfile-argv[2]; 
n-atoi(argv(3]); 
nprime=ato1 (argv[4]); 
l=atəei(arg%x 1); 
reduce=1; 
} 
1f (argc--9) | /* this ais a cross-fam, with output reduction 
cross=1; 
infilel-argv(l 
infile2=argv[2 
outfile=argv [3 
a] 
E 
] 


` ` - ` 


] 
] 
] 
n=atoi (argv[4]) 
nprime-atoi(argv[5]); 
l=atoi (argv[6]); 
reduce=1; 
} 
/* verify that n, nprime and 1l are powers of 2 iti 
1-n$2; 
if (i!=0){ 
printf fam Lata lectora 
primer ("m nn M calling argument n is not a power of 2Nn"); 
exit (); 
} 


i-nprime$2; 


if (10-094 
printe fam... fatal error mae 
Prime £0 ee calling argument nprime is not a power of 2\n"); 
exit (0s 
) 
j2122; 
if (j!=0){ 
Printf("rfam....fatal errer mn), 
DIJDES(CU NS LI S calling argument L is not. a power of 21m) 
exit(); 
} 
Er * 
/* open input file 1, prepare to get the x signal sample data ay 
pt */ 


ifpl = fopen(intilel,"r"); 

tscanft(ifpl, 7.1.51, 5type,stilem), 

/* verify that file n is greater Thanworwerua eo 
if S(filern uam, 


printf£("£am.. atal error \n ); 

Printer weg inputfilel does not contain enough samples\n"); 
fclose(ifpl); 

exit(); 


} 
/* find p — the number of datasets(nprime long)  */ 
p=((n-nprime)/]); 


/* A 
/* a] 

Jt READ IN DATA SAMPLES = 

705 a, 

/* allocate space to read in the sample values */ 

/ dd 
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Sl-(COMPLEX*)cailoc (n,sizeof(COMPLEX)); 
ıf (sl==NULL) { 
Ber ("fam....ınsufficıent space to allocate si\n"); 
exIt(); 
} 
Del) | 
Dow = 
AW O in the real sample values */ 
/* a 
BESEEr-0; i « n; 1++)( 
fscant(1£pl, "sein", £numreal); 
sl[(i].r-numreal; 
zur) 10.0; 
) 


1| 
> 


else { 


/* ea 
/* read in the complex sample values */ 
ie VE 


for (0a i < n; itt) { 
fscanf(ifpl,"$e %e\n",snumreal, snumımag); 
sl[i].r=numreal; 
sl[i].ı=numimag; 

} 

} 


— E 

/* close the input file #1 n 

Le E 

felese (ifpl); 

pt b 

Et this is a cross-fam go get y samples  */ 

ie Ef 

if (cross--1) í 

px ey 
/* open inputfile2, prepare to get the y signal sample data = 
yat 7 


meee = fopen(infile2,"r"); 

Becant(ifp2,"t1 5i",&type,&file n); 

/* verify that file n is greater than or egual to n */ 
Pe (tile n<n) { 


peer ("fam....fatal error\n"); 
A cos ee inputfilel does not contain enough samplesin"); 
Reose atp; 
exit (); 
} 
ee E 
Ar READ IN Y DATA SAMPLES Ev. 
"s d 
/* allocate space to hold the sample values */ 
— Ey 


yl-(COMPLEX*)calloc(n,sizeof(COMPLEX)); 
if (yl==NULL) { 
print fama. ansufficient spacemto allocate yi\n"); 
exit (); 
} 
Ex */ 
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/* read in the complex sample values */ 
forcu-o e acer 
fscanf(ifp2,"$e  $eMn",&numreal,&numimag); 
yl[i].r-numreal; 
yl[i].i=numımag; 


} 


PK */ 
close the inputfile2 S 
e x, 
ftelosoc(ifp25; 

} /* end rf crossei En 

£* f 

R FORM DATASETS Ey 

foe + 

/* allocate space to hold the p-by-nprime datasets  */ 
= = 


s2=(COMPLEX**)calloc (p, sizeof (COMPLEX*)); 
if (s2==NULL) { 
printf("fam....insufficient space to allccacte r s2xXms 
exit(); 
} 
for E) ste) 
s2 [i] = (COMPLEX*) calloc(nprime, sizeof (COMPLEX) ); 
if (s2[i]==NULL) { 


printf ("fam....insufficient space to allocate 32 9m 
exit (0); 
j 
] 
[= i; 
/* form p datasets of nprime samples from the original sample stream  */ 
Yw 2 


for (i=0; i<p; i++) { 
for (j-0; jsnprime; jo 
k=i*l+j; 
SO < EE 
} 
} 


des fy 

/* if this is a cross-fam form y datasets */ 

ya = 

if (cross==1) { 

/* xy 
> 2% 


/* allocate space to hold the p-by-nprime datasets  */ 
> = 
v2= (COMPLEX**)calloc(p,sizeof (COMPLEX*)); 
if (y2==NULL) { 
printf (" fam. insufficient space Gerallecate yom 7, 
exit (); 
) 
for (i0; i<p; i++){ 
y2[i]=(COMPLEX*)calloc (nprime, sizeof (COMPLEX)); 
if (y2[i]==NULL){ 
printf ("fam....insuffieclient space to tallocate yin). 
exit (); 


bb 
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) 


* * y 
E :orm p datasets of nprime samples from the original sample stream  */ 
fom (1=0; i<p; i++) { 
mere (j=0; j«nprime; j++) { 
k=i*1+3; 
ya k []=y1 [kx]; 
) 
) 


EE end if cross-1 ai 

UN E 

= APPLY WINDOW TO DATASETS er 
VES * 


/* allocate space to hold the window multiplicand values */ 
a = / 

mencow= (COMPLEX*) calloc (nprime, sizeof (COMPLEX) ); 

if (window==NULL) í 


EN rt fam....insufficient space to allocate sl\n"); 

exit (); 

} 
Ji» n 
/* calculate the window values for the "nprime" sample wide window  */ 
E * 


B: u-0; i<nprime; i++) { 
NS (Ewop1*i)/(nprime-1); 
zZ=COS (y); 
Ando w[1].1= 0.54 - (0.46*z); 
emdow[i].i- window[i].r; 


) 


> */ 
/* apply the window to all rows of s2 uid 
S */ 


corma; 1<p; i++) { 
moe 0; J<nprime; j++) { 
Wili e-s2[11][J] £*window[j].r; 
si i>s2[1i]) [3] .i*wimdow[)].1; 


ds E. 

EF this is a cross-fam apply window to y2 */ 

pe E 

BM" cross--1) { 

/* */ 


moe (1=0; i<p; itt) { 
for (J—0; j<mprime; j++) { 
Eu c2 [4 } 2 * window [3] 227; 
2 [3] — 1m] i*window[)].l; 
) 
} 


MW end if cross-l = 

po x 

/* FET EACH DATASET ROW E 
/* = 


7* allocate space to hold rows of data for passing to the FFT routine */ 


b? 
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* * "i 


x=(COMPLEX*)callecinprıme, sizesfleolrrer 
if (x==NULL) { 
print£("fam....1nsutálwcient space to allocate <7); 
ext (i, 
} 
E */ 
'* allocate space to hold rows of results from the FFT reutine */ 
P * 
S3=(COMPLEX**) calloc(p,sizect (COMPLEX S 
if (s3==NULL) { 
praintf("fam....insufficxent space tc lalilocategss3 ur 
xi 
} 
for (1-07 >) <p; eit) { 
$3({1]=(COMPLEX*) calloc(nprime, sizect (COMPLEX; 
if (s3{i]==NULL) { 


print£("fam....insuf£iicient space ro allocate sen. 
exe 
} 
} 
/* er 
/* cake an EET OE each row cf 52 x 
2 a 
zt Set values in arguments sent to fft */ 
direction=1; 
norm=1; 
for (1-0, <p; ir*)( 
L* copy row of s2 into complex array E 


for (j=0; j<nprime; j++){ 
x[512521311315 
) 


jus go get FFT performed y 
É£t (x, nprime, direction, nom); 
/* copy x values into a row Of s3 77 


for (9-0; -j<nprime; J++)1 
s3111191=x 01; 
} 
} 
/* */ 
/* if this is a cross-fam take an FFT of each row of y2 */ 
E = 
if (cross==1) 4d 
E A 
p E 
/* allocate space to hold rows of results from the FFT routine */ 
Ax iar 
y 3=(COMPLEX“*) calilcgetp, sizece (COMPLEX I 
if (y3==NULL) { 
printf("fam....insufficient space to allocate 3 m 
exite); 
} 
for (i=0; i<p; i++){ 
y3[1]=(COMPLEX*)calloc(nprime, sizeof (COMPLEX)); 
if (y3[i]--NULL)( 
printf ("fam... insufficient space to allocate yan 


b8 


2 12.21 1992  fam/fam.c Page 8 


exit (); 


1 


f 
} 
ie =0,; i<p; 1++)( 
PX Copy row of y2 into complex array ER 
PO) 0) 7<nprıme; j++) { 
ae 
} 


t go get FFT performed E 
ES nprime,direction,norm); 
/* copy x values into a row of y3 ER 


For (3=0, j«nprime; j++){ 
ala 191 =x (31; 
) 

} 


je end if cross-l x 

y* E 

Js E 

ft DOWNCONVERT EACH ROW = 
/ vind 


y=) allocate space to hold the downconversion multiplicands  */ 
VO y 
ANACONV=(COMPLEX**)calloc (p, sizeof (COMPLEX*)); 
if (dwnconv==NULL) { 
printf ("fam....insufficient space to allocate dwnconv\n"); 
exit (); 
} 
se aO i<p; i++) { 
awnconv [1]=(COMPLEX*) calloc(nprime, sizeof (COMPLEX) ) ; 
if (dwnconv[i]==NULL) { 


Peimer("fam....insufficıent space to allocate dwnconv\n"); 
Exc); 
) 
) 
/* m 
peeeaownconvert each of the transform sequences a 
E */ 
Z = calculate the down conversion multipliers us 
/ Eu 


mainfac-twopi*l/nprime; 
ESO; 1<p; i++) { 
for (j=0; j<nprime; j++) { 
convfac=i*j*mainfac; 
dwenconv[1][3].r=cos (convfac); 
dunconvi eNi 1i= (-1.0)*sin(convfac); 
} 
} 
/* multiply downconversion factor against each frequency value  */ 
for (i=0; i<p; i++) { 
for (j=0; j<nprime; j++) { 
Areas ala Ae amncony (A lr -ss[i]fj).i*dwnconv[i][jl.ai: 
ss fa] LJ] -i583 [Li] c:dwneonmcu (3) ¿1 sS3[11 13]. .1i*awnconv[i] 131.5; 
s3[i][j].rsnumreal; 
) 
) 


if (crosszz1) | 


b" 
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* “y 

* if this is a cross-fam downceonvere v3) 

/* = 

* multıply downconversion factor agaınst each frequency value */ 

For (1-0, sp, It)! 
for (3=0; j<nprime; 3++)( 

numreal-y3[i][Jl.r*dwnconvlil[]Jl.£ = y3 (2) [qe ewnconw En ps 
y3[i][j1.i->y3[1][jJ].r*dwnconv[i]ljJl. i v3 [3] Dod nc ee eee 
ya 1 e-—numeeau: 


Ft 


r /* end if cross=l */ 
"ke e 


e MULTIPLY COLUMNS E 
y * * y 


/* allocate space to hold the correlation multiply (column multiply) results */ 
i = 
a=nprime*p/2; 
s4= (COMPLEX**)calloc (a, sizeof (COMPLEX*)); 
ıf (s4==NULL) { 
printf("fam....insufficient space to allocate sa nm i, 
exit (); 
} 
for (i120; i«a; 1++)1 
s4[i]=(COMPLEX*) calloc(nprime, sizeof (COMPLEX) ); 
if (s4[i]==NULL) { 
printf("fam....insufficient space to allocate s4\n"); 
exit (); 
} 
) 
a E 
/* allocate space to hold columns of s4 for passing to the FFT routine */ 
i UM 
tempfft-(COMPLEX*)calloc(p,sizeof (COMPLEX)); 
if (tempfft==NULL) { 


printf("fam....insufficient space to allocate temptft 4i by ia pp), 
esito); 
) 

S set values in arguments sent to fft  */ 

direction=1; 

norm=1; 

if (cross--0) { /* auto-fam */ 

i n 

/* multiply colummnsvor 's30*7 

a e 


for (a-nprime-1; a» -1; a--) | 
for (b=0; b<nprime; b++) { 
/* multiply the columns of s3 with other columns of s3 i 
for ]s01 <P JN 
tempfft[j].r-(s3[j]lal.r*s3[I3] Ibl z)*(s3 yl ades SMS r 
tempfft[j)].i=(s3[j][a].i*s3[)] (bJ r (sS Sg NS TEES EIE 
} 


/ £ 
y FFT EACH COLUMN * 
y = 

"s go get FFT performed */ 
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Bus tempfft,p,direction,mnorm); 
A= copy tempfft result values back into a row of s4 = 
meee j—(p/4); 35((3*p/4)-1); g**)1 

E (nprrime-a-i)*p/2)-J73-(p/4); 

s4[c] [b]=tempfft [3]; 

} 


EB end for b=0... */ 
E end for az0... */ 
EE end auto-fam  */ 
else { /* cross-fam */ 
/* */ 
SW iply columns of s3 with y3  */ 
f= a 


FO nprime-l; a> -1; a--) Í 
gar (b=0, b<nprime; b++) | 
/*  multiply the columns of s3 with other columns of y3 a 
tor (j-0; J<p; ++) í 
EE [| 9) J| = (SS 1] [e] usss (3) lo) .2ierls3[j) lal=2*y3[3) [b] 1); 
i (s 3[)][ A] rey3[jl[p].r)=(s3[j)[a].r=y3[3j] [D] A); 
) 


/* E 

Jt EET EACH COLUMN Ey 

We = 
= go get FFT performed El 
BEstempfft,p,direction,norm); 
/* copy tempfft result values back into a row of s4 a 


A (0/4); j«((3*p/4)-1); j**) ( 
C« ( (nprime-a-1) *p/2)*j- (p/4); 
s4[c] [b]>tempfft[j]; 

} 


meee ena for b=0... */ 
EE end for a=0... */ 
EE * end cross-fam  */ 
x E 
OUTPUT */ 
ut E 
halfp=p/2; 


if (reduce--0) ( /* do not reduce the amount of output */ 


if (((argc==8) && (argv[7] [1] ==’ a’ )) 11 ((argc==7) && (argv[6] [1] ==’a’))) | 
/* make an ascii output file of the whole spectral plane */ 
/* open the output file and place header information in it */ 

prp = fopen(outfile, "w"); 

max num alf-nprime*p/2; 

max num f-nprime; 

Ur OED, v1 Sin”, max num alf,max num f£); 

nore) crsmaxenum-alf; 35 1 

For (7-0, Jenprime, JEFA 
ELI ODID. 2e. seA^n”,s4[i][j].r,s4[i][J].i); 


E / 
FM end for i=... = 
fclose (ofp); 
|a end I asciì... */ 


else make a plot sxaf no reduced output file */ 
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/* open the output file and place header information in it */ 


} 


ofp = fopen (outtile, "w"); 
/* place header information in the file */ 


daran ype-2, /* real or comple 

max num alf-nprime*p/2; /* num alpha */ 
alpha min-0.0; ¿alpha min 2 

alpha max=1.0; A 

max_num f=nprime; AMAR par, 

£f mwn all— O. 5, ¿e F'min a 

f imax kaS RONS; ¿E E iman kan 


fprintf (ofp, "%i $i $e $e $1 $e sein", data type, max rumi 6 a bol ikma 
alpha max,max num f,f min all,f max all); 
for (i=0; i<nprime; i++){ /* group of lines */ 
for (js50; j€(5/2); j**)l /* line im che group, 
/* place alpha line header information in file */ 
CUED ali (1 *p/72) 44; /* alpha number */ 
num f-nprime-i; [* DU 7 
f min- -0.5-(0.5'1i/nprimce m n 
£ max-  0452(0.2 554mpprmel Tema i 
fprintf(ofp,;" $i $1 se sen curs alí Mun te, tm re ma 


for (k=1; k<nprime; k++) { /* porius on the mer / 
tempi=(p/2* (nprime-1))-((k-1) *p/2) +3; 
tempJ=k; 


numreal=s4[tempi] [tempj] .r; 
numimag=s4 [tempi] [tempJ] .i; 
fprintf(ofp,"%e $eAn",numreal,numimag); 
l 7y endfor k= ma */ 
7 eud O je = 

Aena: LOL los ir 

/* close the output file */ 

ECLOSe (Sep); 

/* end else ... "4 


) /* end if reduce=0 ... */ 


else ( /* reduce the amount of output  */ 


if (((argc==9)&& (argv[7] [l1]==’a’)) || ((argc==8)&& (argv[6] [1]==’a’))) { 


/* 
yx 


make an ascii output file of the whole spectral plane */ 
open the output file and place header information in it */ 
ofp = fopen(outfile, "w"); 
max num_alf=nprime; 
max num f-nprime; 
fprintf(ofp,"£i ti|n",max num alf;max nume); 
for (i=0; i<max num alt, at 
Eor (j=-0, caprine; J+) 
numreal-s4[i*halfp][j].r; 
numimag-s4[i*halfp][j].i; 
bigmag-sqrt (numreal*numreal - numimag*numimag); 
for (k=1; k<halfp; k++) { /* look for largest value */ 
tempreal-s4[i*halfp-k][jl.r; 
tempimag=s4[i*halfp+k] [j] .i; 
tempmag=sqrt (tempreal*tempreal + tempimag*tempimag) ; 
if (tempmag>bigmag) { /* found a larger value */ 
bigmag=tempmag; 
numreal=tempreal; 
numimag-tempimag; 


TE 
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} /* end if tempmag>bigmag... * 
hb :/* end If kamia ey 
c M »e\n",numreal, numımag); 
EE. end IOEg-... */ 
EE /* end for un = 
ie ose (OLD) ; 
EE end if ascii... */ 


NU. make awplot sxaf reduced output file */ 

7* open the output file and place header information in it */ 
Esp = fopen(outfile,"w"); 
EN cce header information in the file */ 


pura type-2; 4 real or complex */ 
max num alf-nprime; /* num alpha */ 
alpha min-0.0; am alpha min */ 
alpha max-1.0; 7* alpha max */ 

max num f-nprime; 7^ max mum */ 
min _all- Ur “i mina 

f max . 0. n /* f max 2) rf 


fprintf(ofp,"$i $i $e $e $i $e $eWn",data type,max num alf,alpha min, 
alpha max,max num f,f min all,f max all); 
for (i=0; i<nprime; i++){ /* group of lines */ 
/* place alpha line header information in file */ 
Gurr £ /* alpha */ 
num f=nprime-i; /* numeros 
£ min= -0.5+(0.5*i/nprime); f= m y 
memax= 0.5 (0.5*i/nprime); is £ mem */ 
FELIN (Op, $1 Giese seXn" ACUEE alee hum 2, o min, E max); 
for (k=i; k<nprime; k++) { /* points on the line */ 
tempi= (halfp* (nprime-1))-((k-1i)*halfp); 
tempj=k; 
numreal-s4[tempi][tempj].r; 
numimag-s4[tempi][tempj].i 
bigmag=sqrt (numreal*numreal + numimag*numimag) ; 
Poe o j++) Yeu look tor the largest value */ 
tempreal-s4[tempitj][tempj].r 
tempimag=s4 [tempi+j] [tempj] .i 
tempmag-sgrt (tempreal*tempreal + tempımag*tempimag); 
ıf (tempmag>bigmag){ /* found a larger value */ 
bigmag=tempmag; 
numreal=tempreal; 
numimag-tempimag; 


E ena it tempmag> 2.2.2 */ 
lee end. For je. i 
EDI nL (OLD, 4e sein", numreal, numimag); 
} /* end for k=. xf 
IN * end for is... E7 


“close the output file  */ 
Eelose(ofp); 
} /* end else binary... E 
) /* end reduce the amount of output */ 


J /* end of main */ 


?3 


APPENDIX C.  SSCA PROGRAM USE 
Correct commands are: 
tor the cross-ssca 


ssca inputfilel inputfile2 outputfile n nprime 1 Oflag 
ssca inputfilel inputfile2 outputfile n nprime 1 Oflag red 


or for the auto-ssca 


ssca inputfilel outputfile n nprime 1 Oflag 
ssca inputfilel outputfile n nprime 1 Oflag red 


where: inputfilel is the file containing the x signal samples 
inputfile2 is the file containing the y signal samples 
outputfile is where the spectrum values will be placed 
n is the number of samples to use from the inputfiles 
and it must be a power of 2 
nprime is the group size of the input datasets and it 
must be a power of 2 
l is the starting offset of subsequent datasets and 
it must be a power of 2 
Oflag indicates the output filetype, ascii or binary 
-asc for an ascii, MATLAB compatible file, 
full cyclic spectral plane 
-plo for a plot_sxaf, SSPI compatible file, 
half cyclic spectral plane 
red reduces the amount of data output 


note I: the first two entries in every input file are 
expected to be the datatype andini datatype = 1 for real 
values, 2 for complex values. n is the number of samples 


contained in the file, one sample per line. 


note 2: error may occur using the reduce option if 
(n-nprime)/nprime is not an integer. 


input file format: 
datatype n 


sample #1 
sample #2 


sample #n 
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GUEDUESASCIT MATLAB format: 
nprime n 


uem C1) 
ousmput (1) (2) 


output (nprime) (n) 


SUepule Sort plot sxaft format: 
datatype (1 for real, 2 for complex) 
number of alphasalphaminalphamax 
number of freqsfreqminfreqmax 
value of alpha #1 


number of fregs in #1freqmin_in_#1freqmax_in_#1 
output for fregmin in sl 


output for fregmax in #1 


value of alpha talphamax 
number of freqs in #alphamaxfreqminfreqmax 
output for freqmin_in_#alphamax 


output for freqmax_in_#alphamax 
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Dee 


*in 
> al 
rin 
rin 
zin 


OYE 


EE -2201992 —Strip/ssca.c Page i 


exude <stdlib.h> 

ende <stdio.h> 

clude <math.h> 

“de "/home3/carter/thesis/SSPI/pam/fft.c" 
clude "/home3/carter/thesis/SSPI/pam/radix.c" 
main (argc,argv) 


argc; 


Ghar *argv[]; 


{ 


COMP LEX 


xy 2453, SZ WNEONYV:; 


REXUPEPEX *tempfft,*sl,*yl,**s2,*window; 
um y. k,type,n,nprime,p,l,r,direction, norm; 
meme, CLOSS, file n,m, reduce, $1zes3, redindex; 
int tempi, tempj, tempk,templ; 
ESutint,; 
at numreal,numımag, convfac,mainfac; 
at bigmag,tempmag,tempreal,tempimag; 
at twopi-6.28318530718; 
at *outreal; 


E 
Elo 
UO 
Eo 
BO 


double z, y; 


cha 


/* 


/* 


meri lel, ~infile2, *outfile; 

MU ifDl,*ifp2,*ofp; 

X passes data to/from fft routine 

S3 holds results of operations after first fft 
sl holds the x sample values from inputfile 1 
yl holds the y sample values from inputfile 2 
s2 holds the channelized x input data from sl 


ay 


check for the correct number of input arguments 
the correct commands are 
Sscaminputtilewourpurftile mn nprime 1 Oflag 
Ssca inputfile outputfile n nprime 1 Oflag reduce 
ssca inputfilel inputfile2 outputfile n nprime 1 Oflag 
ssca inputfilel inputfile2 outputfile n nprime 1 Oflag reduce 


where: 


mote 1: 


note 2: 


inputfilel is where the x signal samples are expected to be 
inputfile2 is where the y signal samples are expected to be 
outputfile is where the spectrum values will be put 
n is the number of input samples to use 
nprime is the group size of input datasets 
l is the starting offset of subsequent datasets 
Oflag indicates the output file type 
Sas fOr an ascii, matlab compatible file 
Po tar ian tasa. plot sxaft compatible file 
reduce is an option to generate a reduced amount of output 


Fn I aresexpected to be on the first line of the inpuerrle/s 
type — 1 for real, 2 for complex, 
cen obser tot Samples following the first line 


error may occur if (n-nprime) is not evenly divisible by nprime 
when using the reduce option 


maintenance history 
"wrastcenoby LCDR Nancy J. Carter, USN NPGS Monterey CA OL "October 1992 


n 
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15 Oct 92....copied from current version of fam.c to älter to ssca alqorieas 
. . . modified to perform ssca per Brown/Gardner/Loomis paper 
20 Oct 92....made outputfile ascii format MATLAB compatible 
22 Oct 92....made an output format SsPr loteo Sa cts 
....24 Nov 92....added option for a reduced amount of output 
n 
aX * Ht 
NC CHECK INPUT VALUES ay 
ee eM 
if ((argc !- 7)&&(argc !- 8)&&(argc !- 9))( 
printf("ssca...£atal errcr'n 
Et DE Mm incorrect number of calling arguments\n"); 
PEN eee correct formats are:\n"); 
PEENE c cm ssca inputfile outputfile n nprime 1 Oflag\n"); 
pruni lec ssca inputfile outputfile n nprıme 1 Oflag reduce\n"); 
¡a saa ssca inputfilel inputfile2 n outputfile nprime 1 Oflag\n"); 
p rm e F Gua i ssca inputfilel inputfile2 n outputfile nprime 1 Oflag reducen"); 
printer where\n"); 
pL DN a inputfilel contains the x signal samplesin"); 
PDELIHESQU A. D inputfile2 contains the y signal samples Wn"); 
pre un s Gu am aa outputfile will contain the resultsin"); 
p puma E Cu Puan: n is the number of input samples to use\n"); 
Pants a EN: nprime is the group size of input datasetsin"); 
grant e ("n ns 1 is the offset of subsequent datasetsin"); 
prune M c Oflag is the output file format n" 55 
Dine ty. ees Oflag = -asc, a MATLAB file is producedìn"); 
PRL ad Oflag = -plo, a plot sxaf file is produced un)» 
ai t Lee reduce is an option for reduced amount of output\n"); 
EXIE(); 


} 

if (argc==7) (| /* this is an auto-ssca, no output reduction  */ 
cross=0; 
infilel-argv[1]; 
outfile=argv[2]; 
n-atoi(argv[3]); 
nprime-atoi(argv[4]); 
l-atoi(argv[5]); 
reduce=0; 

} 

if ((argc-58)&&(argv[7][0]!2'r')) (| /* this is a cross-ssca, no data reduction */ 
cross=1; 
infilel-argv[1] 
infile2-2arqv [2] 
outfile-argv[3] 
n-atoi(argv[4]) 
nprime=atoi (argv[5}); 
l=atoi(argv[6]); 
reduce=0; 


; 
; 
, 
; 


1 


if ((argc-28)&&(argv[7][0]2-'r')) (| /* this is an auto-ssca, with data reduction */ 
cross=0; 
infilel-argv[1]; 
outfile-argv [2]; 
n=atoi largv[3]); 
nprime=atoi(argv[4]); 
l=atoi (argv[5)); 
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reduce=1; 


BENE Grgce--9) | ,* this is a cross-ssca, with data reduction ~ 
cross=1; 
infilel-argv[1]; 
Es3re2-argv [2]; 
outfile-argv[3]; 
n-atoi(argv[4]); 
nprime=ato1(argv[5]); 
l=atoi(argv[6]); 
reduce=1; 
! 


EN crify that n, nprime and L are powers of z i: 
i=n%2; 
fee 2! =O) { 

Damer (C ssca...fatal error\n"); 

EE (o... calling arqument artis not a power of 211); 

exit (); 


! 


i=nprime%2; 


rape =O) í 
mente ("Ssca...fatal error\n"); 
PEL... calling argument nprime is not a power of 21n"); 
exit (); 
} 
j21$2; 
nud 0)1 
EpEuBEf("sscaW*.fatal errorWMn"); 
Eun E(".. ll... callino arqument 1 15 /n0t.a pewes Of 24n"); 
exit (); 
} 
/* a 
/* open input file 1, prepare to get the x signal sample data d 
yr E 


EE £open(infilel,"r"); 
een (itpl, "s1 $1",¿type, ¿file n); 
EN: that file n is at least equal ton */ 
ile n < n)( 
FED Ssca...fatal errorAn"); 
IE o inputfilel does not contain enough samples\n"); 
relese(ifpl); 
exit: (); 
) 
End p — the number of datasets (nprime long)  */ 
pz» ( (nm nprime)/1); 
/* find one dimension of final data array */ 
Sizes3-p*1; 


Aot Ef 
Br READ IN DATA SAMPLES */ 
a E 


/* allocate space to read in the x sample values */ 

m t 

SUS (COMPLEX*)calloc(n,s1zeof (COMPLEX)) ; 

if (s1l==NULL) { 
Prine SSsca.-INsutticient space to allocate si\n"); 
exit (); 


T1553 
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if (type--1) | “* read in the real sample values */ 
for Mi=07) i<nA IT )098 
fscäanf (ifpl, "Sein! snumrea 
sl[i] .rsnumreal; 
sl[r].x-0:0; 
! 
} 
else { /* read in the complex sample values */ 
for (150; 3 «np DE 
fscanf(ifpl,"$e  £eWMn",&numreal,é&numimag); 
s1[i].r=numreal; 
s1[1].i=numimag; 


1 
> 


} 


/* */ 

/* close the input file #1 N 

/* */ 

fclose(ifpl); 

/* x 

if (cross==]) { 

/* t 
/* open inputfile2, prepare to get the wy sen ec um Ue s 
/* * / 


itp2 — fopen(infile2, r5 
fscanf(fp2," $1 $1",8tVDCOS TO eum, 
/* verify that file n is orale asp amp en "n 
Ee eer < mye} 
Printt ('ssca...teataleerron ma 
pw rt F GE na sss: inputfile2 does not contain enough samplesin"); 
pelose (G fp2); 
exit (); 
} 
/* */ 
/* allocate space to read in the y sample coles. 
ie zu 
y1=(COMPLEX*)calloc (n, sizeof (COMPLEX)); 
if (yl==NULL) { 
printf ("ssca...insufticient space tojadlócarte ini 


exit(); 

) 
/* 27 
/* read in the complex sample values */ 
j* x 


for (i-0, i<. mn; 1++)1 
Escant(ifpe ss se\n", &numreal, &numimag) ; 
y1[i].r=numreal; 
y1[i].i=numimag; 

) 


jk */ 
/* close the inputfile2 * / 
/* x 


fclose(ifp2); 
] PZ s end if eGEcSS el, * 
Js n 
PE FORM DATASETS x 
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_* * / 
1 allocate space to hold the p-by-nprime datasets  */ 
NER */ 


ESSICOMPLEX**)calloc(p,sizeof (COMPLEX*)); 
ıf (s2==NULL) { 
prentf("ssca...insufficient space to allocate s2\n"); 
exit); 
} 
BOLE-0. i«p; i-c-t)( 
s2[1]=(COMPLEX*)calloc (nprime, sizeof (COMPLEX)); 
if (s2[i]==NULL) { 
MImenci(’Sssca...insufficient space to allocate s2\n"); 


Exuat (); 
} 
) 
pe */ 
/* form p datasets of nprime samples from the original sample stream  */ 
Mo = 


ton G0; i<p; i++) { 
10-0; j«nprime; j*-*)1! 
ilr); 
s2[i][j]-si(k]; 
} 
} 


es mes 
hee APPLY WINDOW TO DATASETS E 
g5 Bd 


/* allocate space to hold the window multiplicand values */ 
7 E 
window- (COMPLEX*)calloc (nprime,sizeof (COMPLEX)); 
ıf (window==NULL) { 
printf("ssca...insufficient space to allocate s1Mn"); 
exit (); 


} 


i= S 
/* calculate the window values for the "nprime" sample wide window  */ 
ie BA 


Boro1-0; i«nprime; i-*-*)í 
y-(twopi*i)/(nprime-l); 
z-cos(y); 
windew[i].r= 0.54 = (0.46*z); 
window[ti!.r- window[i].r; 


) 


u^. */ 
/* apply the window to rows of s2 containing data D 
Ma En 


for (1=0; i<p; i++) { 
for (j=0; j<nprime; j++) { 
2 |a] [3 J .22Ss2 [si] (3 ] os *window[ 7} .1; 
scm scssDy  r*windowl)j].r; 


jut =, 
/* FFT EACH DATASET ROW a 
s ui 


/* allocate space to hold rows of s2 for passing to the FFT routine */ 
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* x 
x= (COMPLEX*)calloc (nprime, sizeof (COMPLEX)); 
1f (x==NULL) { 
printf("ssca...insufficient space to allocate = mr 
exit (); 
} 
x = 
/* allocate space to hold rows of results from tae reinas 7), 
ix x 
s3=(COMPLEX**)calloc(nprime, s1ze0f (COMPLE DD 
if (s3==NULL) { 
print£("ssca...insufficient Space tS oca Ir 
exit(); 
} 
tor (1-0, Tr_enprime, ir): 
s3[i]-(COMBLEX*)calloc((p* 1 seen WET O'N 
if (s3[i]==NULL) { 
printf("ssca...insufficient space to allocate s3\n"); 


exit (); 

} 
} 
Ax i 
/* take an FFT of each row Of S2 27 
(> T 
A set values in arguments sent to fft  */ 
direction=1; 
norm-1; 
tore A d 

[= copy row of s2 into complex array i£ 


tor (53-05 j<nprime att) 

xs 39, 

) 
Bi go get FFT performed =” 
£ffe (x, nprime,direetion, nern); 
/* Copy X values into Ich column of ss 5 
r=1 51; 
for (j=0; j<nprime; j++){ 

Soi 2 
} 
} 


JE d 
i DOWNCONVERT EACH ROW  */ 
yw * y 


/* allocate space to hold the downconversion multiplicands  */ 
/* */ 
dwnconv-(COMPLEX**)calloc(p,sizeotf (COMPBEEX*)5; 
if (dwnconv==NULL) { 
printf("ssca...insufficient space to allocate dwnconv\n"); 
exit (); 
} 
for (i=0; i<p; 144) 
dwnconv [i] =(COMPLEX*) calloc (nprime, sizeof (COMPLEX) ); 
if (dwnconv[i]--NULL)í 
printf("ssca...insufficient space to allocate dwnconvin"); 
exit (); 
} 


Be 
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downconvert each of the transform sequences * | 
7i 
calculate the down conversion multipliers St 
we 


KH 


Mainfac=twopi*l/nprime; 
for (m=0; m<p; m+t+){ 
por (k=0; k<nprime; k++) { 
convfac=mainfac* (k- (nprime/2)-71)*m; 
dwnconv [m] [k] .r=cos (convfac) ; 
ECXconvim][k].1— (-1.0)*sin(convfac); 
) 
) 
7* multiply downconversion factor against each frequency value  */ 
(ic; 1<p; i++) { 
tempi-1*1; 
meme jJ=O; j<nprime; j++) { 
mumceadl=s3[ 1) [(temp1] .r*dwnconv[i][3].x -— s3[j] [tempi] .1*dwnconv[i] [j).1; 
eee temnp1].1353[3] [temp1] .r*awnconv[1][3].1 + s3[ 3] [(temp1] .1*dwnconv [1] [3].x; 
s3[j] [tempi].r2numreal; 
} 
} 


Er Ev 
or REPLICATE COLUMNS s 
ye E 
Rm 1-0; i«p; itt) Í 
fed. * |; 


pewcopy column into 1-1 adjacent columns */ 
bore (j—0; j<nprime; j++) { 
for (kei; k<l; k++) { 
SET sut], 


INS end for ke... Ey 

EM * end for je...  */ 
EE end for i=... */ 
e n 
en MULTIPLY COLUMNS i 
kt 2 
BE c:055--0) | /* auto-ssca  */ 
Z E 
Zxmultiply columns of s3 with conjugate value of si */ 
i E 


for (a=0; a<sizes3; a++) { 
for (j=0; j«nprime; j++){ 
Se sxbLae*stialsr)-(ts3[3] [a] .i*si[a].i); 
Soha (s'i alu si[al.c)-—(s3[J] [la] .r*sl[a]. 1); 
) 
mer end for a=0...  */ 
) /* end auto-ssca */ 
else ( /* cross-ssca  */ 


(ee E 
"multiply columns of s3 with conjugate value of yl */ 
yt y 


for (a=0; a<sizes3; a++) ( 
ien suprime; Al 
So ai pss (s i aman ym ann; 
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s3[3] [a].32(83[3] [a] -x*yd [a ox) SOS SI Ia EGO ES EE 
) 

Ll, "end. 5r a-O. XE 

/* end cross—ssca “Y 
* */ 
ee FFT EACH ROW x 
/* */ 
* x2 


allocate space to hold rows of s3 for passing tô thew route 
tempffit=(COMPLEX*) calloc(sizes3, Sizeof (COMPLEX) ); 
if (tempfft==NULL) { 
printf ("ssca...insufficient space to allocate tempttt si by A'n, sizes 
exit(); 
} 
i set values in arguments sent to fft */ 
direction=1; 
norm=1; 
for (j=0; j<nprime; 44774 
/* copy row of s3 Into temperte y 
for (k=0; k<sizes3; k++) { 
tempfft ([k]=s3[j] [k]; 
} 


i go get FFT performed ah 
££t (tempfit, sizes3, direction, norm); 
/* copy tempfft result values back into a row of s3 - 


for (k=0; k<sizes3; k++) 
s3[j] [k]=tempfft [k]; 
} 
J Send For =o m F< 


/* */ 
L* ,OUTEUTL xA 
/* */ 


if (reduce==0) { /* no output reduction */ 
if (((argc--7)&& (argv[6] (1]2-'a'))l|| ((arge--26) && (argv[5] [1] 22*a" ) )) { 
/ make an ascii output file full spectral plane 
/* open the output file and place header information in it */ 
ofp = fopen(outfile,"w"); 
fprintf£f(ofp, i SiAn", nprime sizes o 
/* write out data values to file for plotting | 7 
for (i=0; i<nprıme; i++){ 
for (jmo, ISSIZESI) IFA 
ftprintt(otp "$e. Sev’; sata). so ela eee, 
POPE a E 
hr. Se ro.) 
/* close the output file =*/ 
CLOSE (CTD) 
y ALEA IS CA EU 
else ( /* make a plot_sxaf output file half spectral plane */ 
/* open the output file and place header information in it */ 
ofp - fopen(outfile,"w"); 
a=0; /* alpha = 0 */ 
fprintf(ofp,"%i $i te $e $1 $e %e\n",2,51z653,0.0,1.0,nprime,—0.5,9.35); 
for (150; i«nprime; i++) { /* group of lines */ 
for (j=0; j<(sizes3/nprime); j++) ( /* line in the group */ 
tempreal= -0.5+(0.5*i/nprime); 


au 
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tempımag= 0.5-(0.5*ı/nprıme); 
inca (Oi 2051559 Sen ,a, (APrime=1), ,tempreal,tempimag)'; 
a=a+1; 
COKRO; kSkaprime-1); k++) { f* points om the lane */ 
tempi=nprime-k-1; 
tempj-(irtk)*sizes3/nprimetj; 
numreal-s3[tempi] [(tempjJ].r; 
numimag=s3 [temp1] [temp3].1; 
far inef(ofp, "se se\n”,.numreal, numimag); 


nd for ksi k 
m end for gJs... 7 
MW end for 1=... */ 


ye close the output file */ 
fclose (ofp); 
Bir end if else */ 
EE end if reduce-z0... */ 
1f (reduce--1) (| /* reduce the amount of output data */ 
Tf (((argc--8)&& (argv[6] [1]22'a')) || ((arzgc227) && (argv[5] [1]22'a'))) { 
make an ascıi output file full spectral plane */ 
/* open the output file and place header ınformation ın ıt */ 
ofp = fopen (outfile,"w"); 
Rent i(otp, "$1 %i\n”,nprime,p); 
redindex=1; 
Pewee write out data values to file for plotting  */ 
for (i1=0; i<nprime; i++) { 
for (j=0; J<p; jt+){ 
numreal-s3[i][j*redindex].r; 
numimag=s3[i] [j*redindex] .i; 
bigmag=sqrt (numreal*numreal + numimag*numimag) ; 
for (k=1; k<redindex; k++) { /* look for largest value */ 
tempreal=s3[i] [j*redindex +k].r; 
tempimag=s3[1i] [j*redindex +k] .i; 
tempmag=sqrt (tempreal*tempreal + tempımag*tempımag); 
if (tempmag>bigmag) { /* found a bigger value, swap */ 
bigmag-tempmag; 
numreal-tempreal; 
numimag-tempimag; 
EE end it tempmag>.... */ 
|] *cend for kel... */ 
print t(ofp, "se  se\n”, numreal, numimag) ; 
tor je)... t7 
EN er 0D. o t 
E close the output file  */ 
Belose (op); 
EN Ut ascril.»t 


else { /* make a plot sxaf output file with reduction half spectral plane */ 


/* open the output file and place header information in it */ 
ofp = fopen(outfile, "w"); 
redindex=sizes3/nprime; 
a=0; /* alpha = O */ 
rs amis (o E pus Wc Se Je sr Se 2eAYn"',2,nprime,0.0,1.0,nprime,—0.5,0.5); 
for (i=0; i<nprime; i++) [ /* nprime alpha levels */ 
tempreal= -0.5+(0.5*i/nprime); 
Gem-'mag— 0.5—(0.5*i/nprime); 
fprintf(ofp,"$i $i %e %e\n",a, (nprime-i),tempreal, tempimag) ; 
a=a+1; 
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ror (k=0; k< (nprime 1) KRE "* points atWNebhs alpha level 
tempi-nprime-k-1; 
tempj=(1+k) *redindex; 
numreal=s3[temp1] [temp3].xr; 
numimag=s3 [temp1] [temp3].i; 
bigmag=sqrt (numreal*numreal + numımag*numimag); 
for (j=l; j<redindex; Jt) (m line in Er serolp 
tempreal-s3[tempi] [tempjJ +j] r; 
tempimaçg=s3[tempi] [tempj +j] .1; 
tempmag-sqrt (tempreal*tempreal + tempımag*tempimag); 
if (tempmag>bigmag) { /* found a bigger value, swap */ 
bigmag=t empmag; 
numreal=tempreal; 
numimag-tempimag; 


| /* end if tempmad- SS 
H end roro s m 
fprintf(ofp,"3e $e Vn" »nmumreald o naumrmag), 
| /*.end for ke... 09 
k. /* endctoro1m c 77 


/* close the output Ellen 
fclose(ofp); 
) of" cnd irtSelse = / 
|] s end 6 Yr equece= 1 n 
} 
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APPENDIX E. SUBFAM PROGRAM USE 
Correct command is: 


subfam n inputfilel offseti freqzerol conjl inputfile2 offset2 
freqzero2 conj2 outputfile IOflag Itypes 


where: 
n is the number of samples to use from each inputfile, 
it must be a power of 2 


inputfilel is the first file containing signal samples 


offsetl is the initial sample position offset location 
MW inputrilel 


fregzerol is the flag indicating which freguency 
modit icat Iiom to perform on inputfilel, options are: 
zn - zero neg freqs, shift pos freqs down 
zp - zero pos freqs, shift neg freqs up 
zz - don't do anything to frequency components 


com iegenesrlagzsingreating Gf Conjugation or the 
data from inputfilel is to be performed prior to 
multıplicatıon with data from inputtilez, options: 
y - yes, conjugate data 
n - do not conjugate data 


inputfile2 is the second file containing signal 
samples 


offset2 is the initial sample position offset location 
in inputfile2 


freqzero2 is the flag indicating which frequency 
med I a o o pertorm on dinputfile2, options are: 
zn - zero neg freqs, shift pos freqs down 
Zp - zero pos freqs, shift neg freqs up 
ZZ - don't do anything to frequency components 


cen che flag indicating 1f conjugation of the 
data from inputfile2 is to be performed prior to 
cicle clon with data From inputfilel, options: 
y - yes, conjugate data 
n - do not conjugate data 


outputfile is where the spectrum values will be placed 
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IOflag indicates the datatypes of inputfilel, 
inputfile2 and the outputfile. three characters are 
given to indicate the datatype of each file. 
valid datatypes are: 

a - file is of type ascii 

b - file is of type binary 

ex. "abb" means inputfilel is of type ascii, 
inputfile2 is of type binary and the outputfile 
is to be of type binary 


Itypes indicates the types of data inside inputfilel 
and inputfile2. the outputfile is always complex 
floating point, real and imaginary components on the 
same line. two letters must be given to indicate 
the datatype for each input file. all types are 
converted internally to complex floating point. 
valid types are: 

r - Single real floating point samples per line 

i - single integer sample per line 

C - two real floating point numbers per line, 
representing the real and imaginary parts 
of a single sample 

ex. "ic" means inputfilel has integer samples 

and inputfile2 has complex floating point samples 


example: 


subfam 8192 datal.dat 64 zn n data2.dat 2156 zp y 
outputA.dat bbb rr 


input file format: gutput Ea FC OE Ee 
sample #1 magnitude #1 phase #1 
sample #2 magnitude #2 phase #2 
sample #n magnitude #n phase #n 
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aınelude <stalib. h> 

sinclude <stdio.h> 

*include «math.n» 

include "/home3/carter/thesis/SSPI/pam/fft.c" 
sinclude "/home3/carter/thesis/SSPI/pam/radix.c" 
/* values used in up/down conversion */ 


$define INTFREO 49375000.0 /* intermediate signal freq an Hz */ 
sdefine SAMPTIME 0.0000051282 /* orig signal sampling period in sec </ 
main (argc, argv) 

int argo; 


char *arGv[]; 

{ 

COMPHEX “sl, 52 555. 

int i,n,direction,norm,öffl,ofr£2, marterın, scaninreger, reach gc 

float tempreal,tempimag, outputmag, outputphase,numreal,intfreg,t,convfac; 

float *readreal,*readimag, *y; 

char *infilel, *infile2, *outfile, *conjl, *conj2, *freqzerol, *fregzero2, “iorlag, 28 7 75 

EILE trol, 35592. 0; 

Ex 
subfam is a program which uses signal samples from two input files to calculate 
spectral values in a single diamond of the frequency-alpha plane 


the command is 


subfam n inputfilel offsetl fregzerol eonjl inpuerılez orrserzzrregze 
conj2 outputfile IOflag Itypes 


n is the number of samples to use from each inputfile, must be pwr of 2 


inputfilel is the name of the first input file containing signal samples 
offsetl is the initial sample position offset location in inputfilel 
freqzerol is the flag indicating which frequency modification to perform 
on inpubfulelevalid options sere: 
zn — Zero negative frequency components Sn post tos 
components down 
zp — zero positive frequency components, shift negative 
components up 
zz — don't do anything to frequency components 
conJl is the flag indicating if conjugation of the data from inputfiled 
to be performed prior to multiplication with data from ınputfile2. 
valid options are: 
y — yes conjugate data 
n — do not conjugate data 


inputfile2 is the name of the second input file containing signal samples 
offset2 is the initial sample position offset location in inputfile2 
freqzero2 is the flag indicating which frequency modification to perform 
on inputfile2, valid options are: 
zn — zero negative frequency components msi tt posse 
components down 
Zp — zero positive frequency components, shi bEeHcdgdbue 
components up 
ZZ — don't do anything to frequency components 
conj2 is the flag indicating if conjugation of the data from inputfile2 is 
to be performed prior to multiplication with data from inputfilel. 
valid Optra sas ane: 


40 
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y — yes conjugate data 
nodo Mot conjugate data 


outputfile is the filename where the spectrum values will be placed 


Torlag ındreartes the datatypes o£ inputfilel, inputfile2 and outputfile. 
three letters must be given to indicate the datatypes for all files. 
valid datatypes are: 

a le 15 Of type ascii 

b — file is of type binary 

ex. "abb" means :nputfiiel 1s of Gype ascii, 
inpherile2 15 or type binary, and 
the outputfile is to be of type binary 


Itypes indicates the number types of data inside inputfilel and 
inputtilez. the outputfile 15 always complex floating point, 
real and imaginary components on the same line. two letters 
must be given to indicate the number type for each input file. 
all number types are converted to complex floating point. 
valid number types are: 

r — single real floating point sample per line 

1 — single integer sample per line 

C - two floating point numbers per sample on each line 

ex. "ic" means inputfilel has integer samples and 
inputfile2 has complex floating point samples 


sample useage: 
Enprames192.datal.dat 64 zn.n data2.dat 2156 zp y outputA.dat bbb rr 


maintenance history 
created 01 October 1992, LCDR Nancy J. Carter, USN at NPGS Monterey CA 
mOs15.92 — added fregzero options - njc 
mw''c:92.— addedserror check on conjl and conjJ2 - njc 
10.25.92 — changed down/up shifting to same code as in old version 


nem ncm qu uq pou s uc A a a a H */ 
/* check for the correct number of input arguments B 
gee (argo != 13){ 
Prinet( Subfam...fatal error...incorrect number of calling argumentsin"); 
DAMET (ONNI); 
(o. correct format is; An"); 
panes 2sub£fam n 3nputfilel offsetl freqzerol conjl inputfile2 offset2 An"); 
pomer (a Ese zero 2 conj2 outputfile IOflaq Itypes in"); 
Printf("\n"); 
feerntt ("\n"); 
Pac". o... n is the number of samples to use from each inputfile, must be pwr of 21n" 
puuntf(".l.... inputfilel is a file containing the signal samples Wn"); 
A offsetl is the initial sample position offset in inputfilel'n"); 
PE freqzerol indicates type of freq zeroing and shifting to do on inputfilelV 
pruntf(" zn....zero negative freqs and shift down Wn"); 
printet (i ZP Zero positive tregs and shift upin"); 
print£f(" 22 dosnosszero'anvshing, do not shift in"); 
Printer... conjl indicates if conjugation of data from inputfilel is\n"); 
Prine t ( desired before multiplication\n"); 
Primer" y....yes conjugate data\n"); 


q1 
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Printer! n....no do not conjugate data ary 
pre usual setup iS conl =m cen — Qua; 
pram EU od inputfile2 is a file containing the signal samplesìân"); 
OE ees. es offset2 is the initial sample position offset in inputrive2 wae 
DEIN LE freqzero2 indicates type of freq zeroing and shifting to do on inputfilez 
rum te E (u zn....zero negative freqs and shift down\n"); 
PENCE 2P....zero positive freqs andgshi SE Up) 
Printer ZZ2....do not zero anything, demuecweshis&kin' 
prone ue conj2 indicates if conjugation ofwdatawtrom inputfile2 :sXn 55 
prrnpntr( desired before multiplication n"); 
pene Y... .yes conjugate data nm». 
PRINCES n....no do not comgugate data'n'" 5; 
Printer usual setup is ‘congjl =n, conj2 — yu»): 
¡a m outputfile 1s where the spectrum values will be placed\n"); 
print m cr IOflag indicates datatypes of inputfilel,inputfile2 and outputfile\n"); 
PEE AAA m e oe ascii,ascii,ascii\n"); 
print eu aab. eU nn ascii,ascii,binary mo 
printer. aa. —— EN. ascii, Dinary, ascii); 
print.” abb. c e Me mo ME asCli,bindry, binary a); 
prompt Baa uc I w A Ama ascıi,ascii,ascii\n"); 
printed bab. Ec AME I ascii,ascudw, binary mom 
PEINE BRA... sm binary, binary,ascli am), 
Printer” bbb... cc na binary, binary, binary\n") ; 
PEENE ares, oot Itypes indicates number types of input samples\n"); 
printer, Lge ee eee Single real float\n"); 
print f (0 3o ON single integer\n"); 
prime” o OA complex, two floats\n"); 
print E'n) 
printf("....sample useage:\n"); 
printf(" subfam 8192 datal.dat 64 zn n data2.dat 2156 zp y cutputA.dat bbb EF 


exit(); 

} 
nsatoi(argvil])ys 
infilel-argv[2]; 


offl=atoi(argv[3]);/ 


freqzerol-argv [4]; 
conjl-argv[5]; 
ınfile2=argv[6]; 


off2=atoi(argv[7]); 


freqzero2=argv [8]; 
conj2=argv [9]; 
outfile=argv [10]; 
10flag=argyv [EE]; 
itypes=argv [12]; 


/* verify that n is a power of 2 Er 
i-n$2; 
LE MASON 
printf ("subtam a Ea error ms 
peintre ap calling argument n is not a power of 2\n"); 


exit(); 


} 


/* verify that conjl SW ares sora 


zy 


1£ (((conj1[0]1=" y”)£¿£ (con31 [07 ="n 0) 1 ((conj2 10112 Y Econo DE 


print 6 (G Subprfamss 


fatal error mn 


e c conjl and conj2 must be set tComwror a nie 


a= 
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exit(); 
} 
AN mf that itypes are r, i or c E 
if (((itypes[0O]!s»'r')&&(itypes[0O]!-2'i')&& (itypes[0]! 
((itypes[1]!-'r')&&(itypes[1]!-'i')&&(itypes(1]!-' 
Pemer("subfam...fatal error\n"); 


PEE ("....... Toes must De set tor, * Or cAn"); 
exit (); 
} 

/* */ 

mm ocate space for the 1-D data arrays  */ 

[* Lid à 


s1=(COMPLEX*)calloc (n, sizeof (COMPLEX)); 

if (sl==NULL)í 
printf("subfam...insufficient space to allocate sl'*n"); 
exit (); 
} 

SU = COMP LEX*)calloc(n,sizeof(COMPLEX)); 

ıf (s2==NULL) { 
Frinef("subfam.. .insufficient space to allocate s2\n"); 
exit (); 


/* allocate space to read in binary data uU 
mm cat )malloc(l*sizeof(float)); 
p umnrt*)malloc(l*sizeof(int)); 


yo */ 
INPUT */ 
PAX a, 
/* open input files, get the signal samples Es 
a = 
Euch (ioflag[0]) í 
case 'a': (Input ciel is an ASCII file */ 
MIU fpl — fopen(infilel," r")) == NULL) í 
EE subiam..c unable to open ascii ınputfilel\n"); 
exit (); 


} 
/* move to the desired starting offset position in the file */ 
for (1=0; 1<off1;i++) { 
Switch (itypes[0]) { 
case 'r': /* single real float per line */ 
if (fscanf(ifpl, "%e", &tempreal) !=1) { 
printf ("subfam...EOF encountered while attempting to reach offsetl\n"); 
ferose(irpl); 
exit (); 
} 
break; 
case ’i’: /* single integer per line */ 
es came ol, 1" Sscaninteger)!—l)( 
printf("subfam...EOF encountered while attempting to reach offsetl\n"); 
fclose (ifpl); 
exieri), 
} 
break; 
default: /* complex, two floats per line */ 
if (fscanf(ifpl,"%e %e", &tempreal, &tempimag) !=2) { 
printf ("subfam...EOF encountered while attempting to reach offset1\n"); 


aa 
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fclose(ifpl); 
exit Oy, 
} 
|| /* end swateb 3typesi 
} /* end IE 1-0, 1<ertl 7 


/* read in n data samples */ 
for (150;1«nj51**)] 
switch (itypes[0]) { 
case 'r': /* single real float penance, 
if (fscan£f(ifpl, $e" Ss — 
printf("subfam...EOF reached while reading inputfilelNn"); 
fclose(ifpl); 
exit(); 
} 
break; 
case 'i': /* single integer per line */ 
1f (fscanf(itpl, "+1" scaninbeger) IN 
printf ("subfam...EOF reached while reading inputfilel\n"); 
fclose(repi); 
exc 
) 
s1[i].r=scaninteger; 
break; 
default: /* complex, two floats per lines, 
if (fscantf(ifpl,"2e te", ssl lei, ss (eee ey 
printf("subfam...EOF reached while reading inputfilel\n"); 
relose (ifpl); 
exit (); 
} 
y /* end switch itypes[0] */ 
be “ena cor SU sr. 
fclose(ifpl); 


break; 
case 'b': /* inputfilel is a BINARY file 57 
if (( 1fp1 = fopent(infilel,*rb")) == NULL ) í 
printf ("subfam...unable to open binary inputfilel\n"); 
exit (); 


) 
/* move to the desired starting offset position in the file */ 
for (is0;i«offl;i-t-*)l 
Switch (itypes[0]) { 
case 'r': /* single real float per Immer’, 
if (rfread(y,sizeof(*y),l FFI — )4 
printf("subfam...EOF encountered while attempting to reach offsetl\n"); 
fclose(rfpl); 
EXIT 
} 
break; 
case ‘i’: /* single integer per line */ 
if (freadiz,sizeot( 2), LEE I 
printf("subfam...EOF encountered while attempting to reach offsetlin"); 
Ecloese (titolo, 
exito 
) 


qu 
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break; 
deraule. complex, two floats per line y 
sra tresdivosazeort(*y),l,lfpl)ltsi)| 
printf("subfam...EOF encountered wnile attempting to reach offsetl\n"); 
EE CU oDi 
exe), 
} 
CAMS ZE EY) tl, fp) Tel) [ 
printf("subfam...EOF encountered while attempting to reach offsetl\n"); 
cT SC YE liye, 
exire; 
} 
} end switch itypes[0] */ 
) /* end if i=0;i<offl = 


/* read inn data samples */ 
for (1=0;i<n;1i++) { 
Switch (itypes[0]) Í 
case 'r': /* single real float per line */ 
cad Ss IzZeot( y), 1,1£p1)!=1)( 
printf("subfam...EOF reached while reading inputfilelin"); 
totosc rts 
exit(); 
} 
sue r= y 
break; 
case 'i': /* single integer per line */ 
BRE Utread(zysrzeotf(*z),l1,1fpl)!s1)[( 
printf ("subfam...EOF reached while reading inputfilel\n"); 
fe rose Gure Du 
exit (); 
} 
A 
break; 
default:  /* complex, two floats per line */ 
if (fread(y,sizeof(*y),l,ifpl)!-1)( 
printf("subfam...EOF reached while reading inputfilelMn"); 
fclose (ifpi); 
exit(); 
} 
sue r= V; 
Bre troadivosi2eoft(*v),l,ifpl)!sl)[ 
printf ("subfam...EOF reached while reading inputfilel\n"); 
fclose(ifpl) ; 
ex OF 
} 
ST tus 
) /* end switch itypes[0] */ 
jy *Xendetor i-0 i«n  */ 
Belsse(irpl); 


break; 

default: ınvalıdaformat specified for inputfilel */ 
r mis (u Sm arm iaia; ailin format specified for inputfilel\n"); 
exit () ; 


jo cendeswuarenerotlag[0]-*/ 


as 
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switch (1oflag[?T]) i 


case 'a': /* inputfile2 is an ASCeTT Erle, 
if (( ifp2 = £fopen (infile2,"rV)) ANUN 
prıntf("subfam...unable to open ascıı ınputfile2\n"); 
exit (); 


} 
/* move to the desired starting offset position in the file */ 
for (mro l<cEE O, ib) | 
switch (itypes[l]) Í 
case 'r': /* single real float per line */ 
if (fscanf(ifp2,"%e", &tempreal) !=1) { 
printf ("subfam...EOF encountered while attempting to reach offset2\n"); 
fclose(rtp2); 
exit(); 
) 
break; 
case 'i': /* single integer per line */ 
if (fscanft(xtp2," $51 ‚sesanntegeris 1), 
printf("subfam...EOF encountered while attempting to reach offset2\n"); 
£cloeeu BO); 
exit(); 
} 
break; 
default: /* complex, two floats per line */ 
if (fscanf(ifp2,"$e %te",&tempreal, &tempimag) !=2) í 
printf("subfam...EOF encountered while attempting to reach offset2\n"); 
Fclosetrtpo); 
exit(); 
} 
bh /* end switch Ieypesa] ^ 
 . end nf 1-0, Sc 9997 


/* read in n data samples */ 
EOL EO, ESA A) d 
Switch (itypes[1]) ( 
case 'r': /* single real float per line */ 
if (fscan£(ifpê, se" «sl '— Ir 
printf("subfam...EOF reached while reading inputfile2ân"); 
relese(r6p2); 
exit (); 
} 
break; 
case ‘i’: /* single integer per line */ 
if (fscanf(ifp2,"$i",scannreger) !—l) I 
printf("subfam...EOF reached while reading inputfile2ân"); 
fclose(1fp2); 
exit (); 
} 
s2[i].r=scaninteger; 
break; 
default: /* complex, two floats per line */ 
if (£scanf(ifp2,"3e $e”, sl es la AS A 
printf ("subfam...EOF reached while reading inputfile2in"); 
fclose(ifp2); 
EXT O) 


} 


jb 
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IN oudeswitcn itypes[l]) */ 
mm end for x-0 1<n  */ 
bclose(rifp2); 


break; 
case 'b': / Tneucrile?“s a BINARY Zile */ 
meee, 1Ep2 = fopen(infilez,"rb")) == NULL ) { 
printf ("subfam...unable to open binary inputfile2\n"); 
exit (); 


} 
/* move to the desired starting offset position in the file */ 
Poem); 1 <OLEZ; i++) { 
switch (itypes[1]) 1 
case “r”: /* single real float per line */ 
if (fread(readreal,sizeof(*readreal),l,ifp2)!-l)í( 
printf("subfam...EOF encountered while attempting to reach offset2\n"); 
felosse(ifp2); 
exit; 
} 
break; 
case ’i /* single integer per line */ 
if (fread(readinteger, sizeof (*readinteger),1,ifp2) !=1) { 
printf (“subfam...EOF encountered while attempting to reach offset2\n"); 
fclose(ifp2); 
exit (); 
} 
break; 
default: /* complex, two floats per line */ 
if (fread(readreal, sizeof (*readreal),1,ifp2) !=1) { 
printf ("subfam...EOF encountered while attempting to reach offset2\n"); 
fclose (ifp2); 
exit (); 
} 
if (fread(readimag, sizeof(*readimag),1l,ifp2) !=1){ 
printf("subfam...EOF encountered while attempting to reach offset2\n"); 
icbose(ifp2); 
exit (); 
} 
ee nd switch types [1] */ 
| end if i-0;i«off2 er 


fe 
. 


/* read in n data samples */ 
for (i=0;i<n;i++)(Í 
switch (itypes[1]) ( 
case 'r': /* single real float per line */ 
if (fread(readreal,sizeof (*readreal),1,ifp2)!=1)( 
printf ("subfam...EOF reached while reading inputfile2in"); 
tclosefp2); 
exit (); 
} 
s2[i].r-*readreal; 
break; 
case 'i': /* single integer per line */ 
if (fread(readinteger,sizeof(*readinteger),l,ifp2)!-1l)í 
printf("subfam...EOF reached while reading inputfile2ìn"); 
telose(1102)% 


q? 
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exit(); 
) 
scaninteger=*readinteger; 
s2[1].r=scaninteger; 
break; 
default:  /* complex, two £leats per line 
if (fread(readreal, sizeof(*readreal),1,ifp2) !=1) { 
printf ("subfam...EOF reached while reading inputfile2\n"); 
TELOSE DIED). 
EXI, 
} 
s2[i].r=*readreal; 
ıf (fread (readimag, sizeof (*readimag) 7 171162) ee 
printf ("subfam...EOF reached while reading inputfile2\n"); 
telaseilirp2), 
exit () ; 
) 
s2[i].i=*readimag; 
} /* end switch Eyes ie 
b. /* end rc —— xcu =) 
telose (TED 


break; 
default: /* invalid format specified tor saputfile28 
print£("subfam...invalid format specified for anputiiie a. 
exit(); 
/* end switch ioflag[1] */ 
/* */ 


/* ZEROIZE REQUESTED FREQUENCY SIDE AND SHIFT APPROPRIATELY */ 
es Ey 
quarter n-n/4; 
switch (freqzero1[1]) { 
case 'n': /* zero negative and shift down */ 
direction-1; /* request a forward transform time-to-freq */ 
norm=0; 
fft (sl, n, direction, norm); /* results returned a sla 
for (1=0; 1<(0/2); 1++)1 /* zero lower half */ 
sl[i].r=0; 
sl[i].i=0; 
} 
direction--1; /* request an inverse transform freq-to-time  */ 
norm=2; 
fft(sl,n,direction,notm); /* xesults returned 4n sc! 
“= downconvert */ 
teg /* starting time vf 
for (1=0 Smp EFE | 
convfac-TWOPI* INTFREQO*t; 
tempreal-cos(convfac); 
tempimag-(-l)*sin(convfac); 
numreal-sl[i].r*tempreal - sl[i].i*tempimag; 
sl[i].i-sL[1].r’tempimas+sl [1] . 1 <eempreal, 
sl[i].r=numreal; 
t=t+SAMPTIME; 
} 
break; 
case 'p': /* zero positive and shift up n 
direction-1;  /* request a forward transform time-to-freq */ 


48 
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norm=0; 
MANSA An recte ron, Norm);  /* results returned in sl  */ 
Bru -H/2)5 i<n; i++)| /* zero upper half */ 
2-0; 
Saati. L=0; 
) 
@irection=-1;  /* request an inverse transform freq-to-time  */ 
norm-2; 
Eu usi nj?;drrection,norm); /* results returned in sl */ 
Exubconvert */ 
t=0.0; /* starting time xu 
p -5(r-—0; i«n; it) { 
convfac-TWOPI*INTFREQ*t; 
tempreal=cos (convfac) ; 
tempimag=sin (convfac) ; 
numreal-si[i].r*tempreal - si[i].i*tempimag; 
si[i] .i=si[i].r*tempimagts1[1i].1*tempreal; 
sl[i].r=numreal; 
Bet +SAMETIME; 
} 


break; 

case 'z': (tao noch una =’ 
break; 

default: /* invalid format specified for fregzerol */ 
printf ("subfam...invalid format specified for fregzerolin"); 
exit (); 


fer end switch */ 
switch (freqzero2[1]) { 


case 'n': /* zero negative and shift down */ 
direction=1; /* request a forward transform time-to-freq */ 
norm=0; 


BEBO? n.direction,norm); "/* results returned in s2  */ 
E30; i<(n/2); it+)! /* zero lower half */ 
s2[i].r=0; 
CUADROS TO 
} 
direction--1;  /* request an inverse transform freq-to-time  */ 
norm-2; 
BEN n.drirection,norm); /* results returned in s2 */ 
Idownconvert */ 
t-0.0; /* starting time UM 
for (1205; i<n; itt) { 
convfac-TWOPI*INTFREQ*t; 
tempreal=cos (convfac) ; 
tempimag=(-1) *sin(convfac) ; 
numreal=s2[i].r*tempreal - s2[i].i*tempimag; 
s2[1].i=s2[i].r*tempimag+s2[i].i*tempreal; 
s2[i].r-numreai; 
t=t+SAMPTIME; 
) 


break; 

case 'p': /* zero positive and shift up */ 
direction-1; /* request a forward transform time-to-freq */ 
norm=0; 


[omn “nn, drrecrion, norm); /* results returned in s2 */ 
for > in i++) 1 /* zero upper half */ 


als 
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Sis 
Ses sumen 


directaion=_1s /* request an inverse transform freq-to-tıme *y 
norm-2; 
fft(s2,n,direction,norm); /* esunbsssreturnocd 3535s v 
"XXDCORhVent et 
t=0.0; /* starting time a 
tor (10s? <m ir I 
convfac-TWOPI*INTFREQ*t; 
tempreal-cos(convfac); 
tempimag-sin(convfac); 
numreal-s2[i].r*tempreal — s2[i] .i*tempimag; 
s2[i].i>s2[i].r*tempimagrs2[i].i*tempreal; 
s2[i].r=numreal; 
t=t+SAMPTIME; 
} 


break; 
case “ze loeo mice 
break; 
default: /* invalid format specified for freqzero2 */ 
printf("subfam...invalid format specified for freqzero2Mn"); 
exit (); 
} /* end switch */ 
yr =/ 
/* MULTIPLICATION  */ 
ja E 
/* create space for the output array Xd 
p “7 


$3=(COMPLEX*) calloc(n, sazeos (GOMBLEEX))) 
if (s3==NULL) { 
printf("subfam...insufficient space to allocate s3Nn"); 
exit (); 
) 
if ((conjl[0]2-'y')&&(conj2[0]2-'y^)) { /* multiply conjugate sl times conjugates zer 
for. (ie0;1*9n;3 4) [í 
sti] r—sll[i]).r*S2[t]. rs so m mm 
s3(1].1=(+15* (s1 61.182 2), yps is m) se eee 
} 
} 
if ((conjl[0]==’y’)é&&(conj2[O0]==’n")) { /* multiply conjugate al times 2 
for O, 0 mn dd) 
s3[1)22=sl la). 27 s2 (1) oars Sa 
s3l1].J—>Sl|]. r*S2[1] 4—sLll' 92 E 
} 
} 
if ((conj1[0]=="n")£6(conj2[0]=="y*)) ( /* multiply sl times eonjügare 7 
for (ISO ESA 
s3a[|i]jur=sumrna mm s> [a mars 0 Ti Es rm pam 
s3 (1) -i=sL (2) 2*s2 (2) eesti ese ie 
} 
} 
if ((conjl[0]2-2'n')&&(conj2[0]9-^n') Y 01-7 * multiply sl @imesss2s 
for (i=0;i<n;i++)( 
s3[i)].rssl[i].r*s2[i] rss do ee, 
s3[i].ı=sl[:i)] ‚2’32[1] ss D] mS DNE, 
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} 
} 
Brake an EFT of n-point s3 Tut 


t; 


it 


= set values in arguments sent to fft  */ 

direction-1;  /* request a forward transform time-to-freq * 
norm-0; 

SaaS n,direction,norm); "7* results returned in s3  */ 

(ies E 

ESUTPUT */ 

a = 


/* take magnitude and phase of s3 and place in the output file * 
EE 
mueca (10flag[2]) { 


case 'a': /* outputfile is an ascayx file */ 
if (( ofp = fopen(outfile,"w")) == NULL ) { 
printf ("subfam...unable to open ascii outputfile\n"); 
exit (); 
} 
/* write out data to the ascıi file = 


for (i=0;i<n; itt) { 
Sep MaG—seet (Sst) .e ss (1)2.5)r(s3lil.i*s3[i) .i)); 
outputphase=atan(s3[i].i/s3[i].r); 
BEGIN (Op, +. I6.6e 16.081”, outputmag, outputphase); 
} 

I close the output file */ 

close (ofp) ; 


break; 
case 'b': 7 Sutpütfile is a binary file */ 

Bell ofp — fopen(outfile,"wb")) 22 NULL ) { 
pnintf('subfam...unable to open binary outputfileWn"); 
exit(); 

) 
/* write out data to the binary file or 


Epor (1=0;1<n;1++) 1 
numEsumgesart((ssSp]-rtss3[Ii]:r)*(s3[i1].i1*s3[i].3)); 
outputphase=atan (s3[i].i/s3[i].r); 

*y=outputmaq; 
me (Ewcite(y,sizeof(*y),1,o0£p)!=1) ( 
Beneerl2subfam...error while writing outputfile\n"); 
telease(örp), 
exit O; 
) 
*y-outputphase; 
FE (were s1zeot (y), 1,0£p)!=1)1 
printf ("subfam...error while writing outputfilein”); 
fclose(ofp); 


exit (); 
} 
NEA end ir i=0 a... */ 
fclose (ofp); 
break; 
default: /* invalid format specified for outputfile */ 
SL IU ub sam... imvalid format specified for outputfile\n"); 
exit(); 


la xcndEswitcheoflag[i2] */ 
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} I end of marna 





BBSK 
FAM 
EE 
ESM 
PAM 
QAM 
QPSK 
SCH 
SSCA 
SSBT 
SUBFAM 


LIST OF ABBREVIATIONS 


binary phase shift keying 

FFT accumulation method 

fast fourier transform 

frequency smoothing method 

pulse amplitude modulation 
quadrature amplitude modulation 
quadrature phase shift keying 
Spectra lmcemrelation function 
Strip spectral correlation algorithm 
Statistical Signal Processing, Inc. 
sub-FFT accumulation method 
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